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ABSTRACT 

Challenges  to  computational  aerothermodynamic  ( CA)  simulation  and  validation  of  hypersonic  flow  over  plan¬ 
etary  entry  vehicles  are  discussed.  Entry,  descent,  and  landing  (EDL)  of  high  mass  to  Mars  is  a  significant 
driver  of  new  simulation  requirements.  These  requirements  include  simulation  of  large  deployable,  flexible 
structures  and  interactions  with  reaction  control  system  (RCS)  and  retro-thruster  jets.  Simulation  of  radiation 
and  ablation  coupled  to  the  flow  solver  continues  to  be  a  high  priority  for  planetary  entry  analyses,  especially 
for  return  to  Earth  and  outer  planet  missions.  Three  research  areas  addressing  these  challenges  are  empha¬ 
sized.  The  first  addresses  the  need  to  obtain  accurate  heating  on  unstructured  tetrahedral  grid  systems  to  take 
advantage  of  flexibility  in  grid  generation  and  grid  adaptation.  A  multi-dimensional  inviscidflux  reconstruction 
algorithm  is  defined  that  is  oriented  with  local  flow  topology  as  opposed  to  grid.  The  second  addresses  coupling 
of  radiation  and  ablation  to  the  hypersonic  flow  solver  — flight-  and  ground-based  data  are  used  to  provide 
limited  validation  of  these  multi-physcis  simulations.  The  third  addresses  the  challenges  of  retro-propulsion 
simulation  and  the  criticality  of  grid  adaptation  in  this  application.  The  evolution  of  C A  to  become  a  tool  for 
innovation  of  EDL  systems  requires  a  successful  resolution  of  these  challenges. 

1.0  INTRODUCTION 

The  ability  to  provide  simulations  of  hypersonic  flows  over  planetary  entry  vehicles  with  well-defined  uncer¬ 
tainties  on  the  predicted  aerothermal  environments  is  critical  to  the  success  of  solar  system  exploration.  As 
mission  needs  evolve  the  demands  on  simulation  of  hypersonic  flow  over  planetary  entry  vehicles  become  more 
complex  and  the  validation  of  those  simulations  becomes  more  challenging.  Consider  the  role  of  computational 
aerothermodynamic s  (CA)  in  the  definition  of  aerothermal  environments  for  representative  entry  vehicles  in 
recent  years[l-3]  and  the  role  it  must  assume  in  the  near  future. 

The  Mars  Pathfinder  program  (entry  1996)  was  the  first  to  use  numerical  simulation  as  the  primary  source 
for  defining  the  aerothermal  environment  for  entry  to  Mars  [4].  The  vehicle  forebody  geometry  was  a  scaled 
Viking  shape.  Viking  entered  the  Martian  atmosphere  from  orbit  with  an  offset  center  of  gravity  and  used  a 
reaction  control  system  to  trim  at  1 1  degrees  angle  of  attack.  In  contrast,  Pathfinder  used  a  spin  stabilized, 
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direct  entry  at  zero  degrees  angle  of  attack.  The  Pathfinder  simulations  were  validated  with  ballistic  range 
data  and  wind-tunnel  data  obtained  for  Viking.  The  pre-flight  discovery  of  a  static  instability  associated  with 
sonic  line  movement  (a  condition  never  identified  in  the  Viking  experimental  database)  provided  a  serendipitous 
opportunity  to  validate  the  prediction  in  flight  through  accelerometer  data.  The  aeroshell  was  not  instrumented 
but  the  accelerometers  recorded  vehicle  oscillations  as  the  sonic  line  jumped  from  the  corner  to  the  nose  and 
then  later  again  as  the  sonic  line  jumped  back  from  the  nose  to  the  shoulder.  These  jumps,  associated  with 
evolving  gas  chemistry,  were  clearly  defined  within  the  matrix  of  simulations  performed  to  define  the  aerodata 
book. 

Support  of  aerothermal  design  for  the  Orion  entry  capsule  and  more  recent  robotic  entry  vehicles  to  Mars 
have  brought  new  simulation  and  validation  needs  to  the  forefront.  Reaction  control  system  jets  fired  from 
the  base  of  the  vehicles  for  active  control  authority  interact  with  the  high  speed  flow  expanding  around  the 
corner  of  the  vehicle  into  the  wake.  [5]  Simulations  are  executed  to  quantify  the  magnitude  of  jet  interaction 
effects  on  aerodynamic  coefficients  and  to  quantify  magnitude  of  any  hot  spots  on  the  aft  heat  shield  that 
may  result  from  these  interactions.  At  present,  such  simulations  are  not  validated.  These  jet  interactions  are 
typically  unsteady.  The  ability  to  adapt  grid  to  the  plume  shocks  and  free  shear  layers  is  limited  by  the  use  of 
structured  grid  topologies.  The  simulations  are  valuable  in  that  model  perturbations  can  be  executed  to  provide  a 
better  understanding  of  the  uncertainties  introduced  by  the  interactions.  But  clearly,  more  advances  are  needed 
if  CA  simulations  are  to  independently  address  these  challenges  with  well  defined  prediction  uncertainties. 
Another  issue  in  current  design  efforts  is  the  prediction  of  transitional  and  turbulent  interactions  associated 
with  tension  ties  and  compression  pads  in  the  Orion  heat  shield.  [6,  7]  These  perturbations  to  a  smooth  thermal 
protection  system  (TPS)  introduce  a  transitional  or  turbulent  wedge  of  flow  downstream.  Transition  criteria 
are  defined  from  ground-based  tests[8]  with  some  corroboration  from  a  shuttle  flight  experiment. [9]  Algebraic 
turbulence  models  still  do  the  best  at  predicting  heating  levels  for  attached,  fully  developed  turbulent  flow  but 
these  methods  are  not  suited  for  simulating  the  transitional  /  turbulent  wedge  behind  such  perturbations. 

Aerothermal  design  support  for  future  vehicles  must  address  new  challenges  introduced  by:  (1)  entry,  de¬ 
scent  and  landing  (EDL)  for  heavy  mass  to  Mars  (one  metric  ton  is  the  approximate  limit  for  existing  archi¬ 
tectures);  and  (2)  limitations  on  payload  mass  to  the  outer  planets  or  return  from  Mars  derived  from  extreme 
heating  levels.  As  will  be  discussed  in  the  next  section  both  of  these  challenges  are  related  to  a  desire  to  mini¬ 
mize  system  mass  by  utilizing  aerodynamic  drag  to  slow  the  vehicle.  In  both  cases,  packaging  more  mass  into 
vehicle  diameters  constrained  by  present  launch  systems  tends  to  increase  the  ballistic  coefficient  of  the  entry 
vehicle.  In  the  case  of  large  mass  to  Mars  if  the  ballistic  coefficient  is  too  large  the  entry  vehicle  runs  out  of 
altitude  before  it  runs  out  of  velocity.  In  the  case  of  high  entry  velocity  (to  outer  planets  or  return  to  Earth) 
a  high  ballistic  coefficient  tends  to  increase  peak  heating  rates  and  requisite  thermal  protection  system  mass. 
The  ballistic  coefficient,  f3  =  m/iCoA)  where  m  is  mass,  Cd  is  drag  coefficient  and  A  is  reference  area, 
can  be  decreased  by  increasing  the  reference  area  with  deployable  structure  (semi-rigid  or  inflatable).  Several 
scenarios  for  deployment  in  the  hypersonic  and  supersonic  domain  have  been  suggested.  Simulation  capability 
must  evolve  to  handle  both  small  and  large  scale  (possibly  unsteady)  structural  deformations  in  defining  the 
aerothermal  environment.  Another  option  is  to  carry  some  retro-propulsion  capability  but  search  for  optimized 
configurations  that  minimize  propulsion  system  mass  (and  possibly  thermal  protection  system  mass)  through 
engagement  in  the  supersonic  (and  possibly  hypersonic)  domains. 

In  the  view  forward  it  is  clear  that  CA  must  be  able  to  simulate  a  variety  of  flow  physics  phenomena  coupled 
with  structural  and  material  response  to  the  aerothermal  environment.  Furthermore,  the  simulation  setup  should 
be  sufficiently  simple  and  robust  to  accommodate  adjoint-driven  optimization  studies  so  that  EDL  architecture 
design  trades  may  be  implemented  using  best-available,  high-fidelity  analyses.  The  simulations  must  resolve 
jet  interactions  with  flow  expanding  into  the  wake  (RCS  interactions)  and  opposing  flow  on  the  windward 
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side  (retro-propulsion).  They  must  accommodate  large  scale  surface  deformation  (inflatable  structures).  These 
topologically  complex  interactions  will  require  robust  grid  adaptation.  The  simulations  must  also  accommodate 
unsteady  flows  associated  with  these  interactions.  We  frame  these  challenges  in  the  context  of  simulation  for 
innovation.  Simulation  for  innovation  has  three  defining  features:  (1)  conceptual  systems  can  be  quickly  set 
up  for  simulation  and  optimization,  (2)  robust  grid  adaptation  to  a  user  specified  discretization  tolerance  is 
available,  and  (3)  credible  estimates  of  simulated  aerothermal  environment  uncertainties  are  published.  Feature 
(1)  is  a  requirement  shared  across  all  speed  domains  but  features  (2)  and  (3)  have  challenges  unique  to  the 
hypersonic  flight  domain.  In  the  material  that  follows  we  will  review  recent  efforts  to  address  the  numerical 
and  physical  modeling  challenges  required  for  CA  simulation  for  innovation  of  planetary  entry  vehicles  and  we 
will  touch  base  with  existing  flight  and  ground-based  data  to  validate  these  simulations. 

2.0  PROBLEM  DEFINITION 

The  suite  of  physical  models  required  to  define  the  aerothermal  environment  of  an  entry  vehicle  are  a  function 
of  relative  velocity  at  entry  interface,  vehicle  size  and  ballistic  coefficient  0  =  m/(C^A)  (A  cannonball 
has  a  very  high  ballistic  coefficient,  a  beach  ball  has  a  very  low  ballistic  coefficient.)  In  general,  planetary 
entry  vehicle  simulations  must  include  capability  to  model  gas-surface  interactions  (catalysis  of  atomic  species 
to  molecular  species),  ablation  (thermal  protection  system  (TPS)  mass  emanating  from  the  boundary  due  to 
material  degradation  at  high  temperatures),  and  radiation  (energy  flux  originating  from  the  high  temperature 
shock  layer).  If  the  vehicle  is  sufficiently  large  the  aerothermal  environment  at  peak  heating  may  be  augmented 
by  a  turbulent  boundary  layer.  In  addition  to  these  processes  long  associated  with  the  simulation  of  aerothermal 
environments  the  demands  of  entry,  descent,  and  landing  at  Mars  have  introduced  new  challenges.  (At  least 
new  to  EDL  -  the  challenges  discussed  below  are  also  relevant  to  innovative  concepts  for  aerocapture.[10])  In 
this  chapter  we  review  the  driving  issues  that  define  the  additional  suite  of  simulation  capabilities  required  for 
computational  aerothermodynamic  (CA)  analyses  of  planetary  entry  vehicles. 

2.1  Entry,  Descent,  and  Landing  —  EDL 

Braun  and  Manning[l  1]  present  EDL  challenges  at  Mars  and  the  need  to  expand  beyond  heritage  systems.  They 
note  that  there  have  been  five  successful  landings  on  Mars  including  Vikings  I  and  II  (1976),  Mars  Pathfinder 
(1997),  and  Mars  Exploration  Rovers  Spirit  and  Opportunity  (2004).  All  five  of  these  successful  systems 
have  landed  masses  of  less  than  0.6  metric  tons.  They  have  landed  at  low  elevation  sites  (below  1  km  Mars 
Orbiter  Laser  Altimeter  (MOLA)).  All  accepted  a  relatively  large  uncertainty  in  landing  location  of  100s  of 
kilometers.  All  of  the  current  Mars  missions  have  relied  on  large  technology  investments  made  in  the  late 
1960s  and  early  1970s  as  part  of  the  Viking  Program.  They  utilize  the  aerodynamic  characterization  of  70- 
deg  sphere  cone,  a  forebody  heatshield  thermal  protection  system  of  SLA-561,  a  supersonic  disk-gap-band 
parachute,  and  autonomous  terminal  descent  propulsion.  These  Viking  heritage  EDL  systems,  the  thin  Martian 
atmosphere,  and  small  scale  height  of  obstacles  on  the  ground  limit  accessible  landing  sites  to  those  below  - 
1.0km  MOLA.  So  far  the  southern  hemisphere  has  been  largely  out  of  reach.  (Approximately  50%  of  the  planet 
surface  remains  inaccessible  with  current  EDL  technologies.)  The  highest  landing  to  date  is  MER-Opportunity 
at  Meridiani  Planum  (-1km  MOLA).  Mars  Science  Lab  (MSL)  is  attempting  to  develop  an  EDL  system  capable 
of  delivering  0.8MT  to  +2km  MOLA. 

The  Mars  Exploration  Program  is  attempting  to  develop  systems  that  deliver  0.8  metric  tons  to  the  surface 
for  MSL  and  1.2  MT  for  Mars  Sample  Return.  In  order  to  prepare  for  human  scale  missions,  we  will  need  to 
develop  EDL  systems  that  can  get  30-100  metric  tons  down  to  the  surface  per  landing.  Viking  heritage  systems 


RTO-EN-AVT-186 


11  -3 


Challenges  to  Computational  Aerothermodynamic 
Simulation  and  Validation  for  Planetary  Entry  Vehicle  Analysis 


ORGANIZATION 


Alt.(km  MO  LA) 


Vel.  (m/s) 


Figure  1:  Phase  plot  for  robotic  MSL  entry  with  ballistic  coefficient  (3 
m/(CDA)  =  100kg/m2  and  lift  to  drag  ratio  L/D  =  0  18.[1 1] 


are  unable  to  achieve  this  goal.  The  problem  can  be  illustrated  using  an  altitude  versus  velocity  phase  chart. 
Figure  1  shows  the  phase  chart  for  the 
robotic  MSL  mission.  Entry  into  the 
sensible  atmosphere  occurs  at  the  right 
boundary.  The  goal  is  to  attain  zero  ve¬ 
locity  at  the  target  landing  altitude  at 
the  lower  left  corner.  The  track  end¬ 
ing  in  a  dashed  line  indicates  the  trajec¬ 
tory  without  supersonic  parachute  de¬ 
ployment  that  hits  the  ground  (0  MOLA) 
at  200  m/s.  The  red  triangular  re¬ 
gion  indicates  the  safe  zone  for  super¬ 
sonic  parachute  deployment  -  bounded 
by  high  and  low  Mach  number,  high  and 
low  dynamic  pressure,  and  minimum 
altitude  requirements.  Viking  heritage 
EDL  requires  a  trajectory  that  enters  this 
triangle,  preferably  near  the  center  of  the 
high  Mach  boundary,  to  enable  super¬ 
sonic  parachute  deployment.  The  track 
ending  in  a  solid  line  in  Fig.  1  indicates 
an  MSL  trajectory  with  supersonic  parachute  deployment. 

Figure  2  focuses  on  the  lower  left  corner  of  Fig.  1  overlayed  with  ballistic  trajectories  as  a  function 
of  ballistic  coefficient  /3.  Increasing  ballistic  coefficient  causes  the  vehicle  to  penetrate  deeper  into  the  at¬ 
mosphere  (lower  altitude)  for  a  given  velocity.  If  the  area  of  the  entry  vehicle  is  constrained  by  launch 
vehicle  diameter  from  Earth  then  ballistic  coefficient  will  increase  with  increasing  payload  mass  to  Mars. 
The  message  from  Fig.  2  is  that  a  ballis¬ 
tic  coefficient  of  150  is  a  limiting  value 
for  current  technology.  Some  relief  is 
attained  with  lifting  bodies  but  the  fun¬ 
damental  limitation  remains  -  heritage 
EDL  systems  cannot  deliver  large  pay- 
loads  to  Mars  with  given  constraints 
of  launch  vehicle  diameter.  Two  new 
EDL  system  elements  are  currently  be¬ 
ing  investigated  to  overcome  this  limita¬ 
tion.  The  first  element  lowers  the  bal¬ 
listic  coefficient  by  deploying  structure 
to  increase  area.  Inflatable  ballutes  de¬ 
ployed  in  the  supersonic  regime  or  prior 
to  entry  interface  are  primary  examples 
of  this  approach.  The  second  element 
augments  aerodynamic  drag  with  retro- 
propulsion  engaged  in  the  supersonic 
and  possibly  hypersonic  regimes. 


Alt.(km  MOLA) 


Varying  /?  from  25-200  kg/m2 
(chute-less  trajectories) 


Figure  2:  Detail  of  Figure  1  showing  influence  of  ballistic  coefficient  on  effec¬ 
tiveness  of  heritage  EDL  systems.  When  j3  is  too  large  altitude  is  too  low  at 
Mach  2  to  deploy  parachute.[11] 


Vel.  (m/s) 
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2.1.1  Supersonic  Retro-Propulsion  —  SRP 

Korzun  et.  al.  [12]  surveyed  experimental  and  computational  efforts  addressing  use  of  retropropulsion  ap¬ 
plicable  to  EDL  systems.  A  single  centered  jet  (emanating  from  the  vehicle  axis)  and  peripheral  jets  were 
compared.  They  report  that  for  low  thrust  coefficients  (ratio  of  jet  thrust  to  freestream  dynamic  pressure  times 
area)  the  peripheral  jet  system  was  superior  to  the  center  in  terms  of  total  axial  force  coefficient  (aerodynamic 
drag  plus  thrust).  The  center  jet  tends  to  create  an  aerospike  effect  that  reduces  aerodynamic  drag.  At  large 
thrust  coefficients  the  aerodynamic  drag  is  overwhelmed  by  the  thrust  and  there  is  no  advantage  to  peripheral  jet 
placement.  As  computational  analyses  mature  with  objective  function  optimization  the  opportunity  to  minimize 
retropropulsion  system  mass  through  tailoring  of  nozzle  design  and  vectoring  should  be  realizable.  Indeed,  en¬ 
gaging  low  thrust  coefficient  retropropulsion  even  in  the  hypersonic  domain  may  be  useful  in  reducing  surface 
heating  if  a  film  cooling  system  can  be  adapted  from  propulsion  system  infrastructure.  The  creative  use  of  CFD 
for  design  in  this  application  is  currently  stymied  by  an  inability  to  effectively  adapt  grid  in  the  jet  interaction 
region  -  even  in  a  supersonic  domain.  An  example  of  this  challenge  will  be  presented  in  the  Numerical  Results. 

2.1.2  IRDT  /  Ballutes 

A  ballute  (balloon-parachute)  or  Inflatable  ReEntry  and  Descent  Technology  (IRDT) [13]  is  an  inflatable,  aero¬ 
dynamic  drag  device  for  lowering  the  ballistic  coefficient  of  planetary  entry  vehicles.  Ballutes  may  be  directly 
attached  to  a  vehicle,  increasing  its  cross-sectional  area  upon  inflation,  or  towed  behind  the  vehicle  as  a  semi¬ 
independent  device  that  can  be  quickly  cut  free  when  the  requisite  change  in  velocity  is  achieved.  Rohrschneider 
et.  al.  [14]  survey  ballute  technology  for 
aerocapture. 

An  example  of  a  trailing  ballute[15,  16] 
is  shown  in  Fig.  3  including  the  computa¬ 
tional  grid  and  color  contours  of  pressure  in 
the  quarter  domain  symmetry  planes.  (Con¬ 
ditions  of  the  simulation  are  for  a  Titan 
Organics  Explorer  with  velocity  equal  to 
8.55  km/s  and  density  equal  to  1.9  x  10-7 
kg/m3. [17,  18]  The  gas  model  includes 
molecular  nitrogen  and  atomic  nitrogen  in 
thermochemical  nonequilibrium.)  An  in¬ 
flatable  tether  connects  the  towing  space¬ 
craft  (lower  right  corner)  to  the  inflated 
toroid.  The  bow  shock  of  the  spacecraft  is 
enveloped  within  the  hole  of  the  toroid,  in¬ 
teracting  with  the  toroid  bow  shock.  A  high 
pressure  zone  in  core  of  the  toroid  drives 
flow  back  toward  the  base  of  the  towing 
spacecraft. 

Ballute  simulation  provides  challenges 
on  many  fronts.  Inflatable  surfaces  deflect 
and  stretch  under  aerodynamic  load  and  rel¬ 
atively  high  temperature.  [19,  20]  Unsteady  interactions  may  follow  from  the  surface  deflections  or  from  a  purely 
fluid  dynamic  origin  (discussed  later  in  more  detail)  and  the  ability  to  accurately  capture  these  phenomena  are 


Figure  3:  Unstructured  grid  used  in  simulation  of  spacecraft  -  tether  -  bal¬ 
lute  system.  Colors  correspond  to  pressure  levels  nondimensionalized  by 
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a  strong  function  of  ability  to  resolve  free  shear  layers  and  the  triple  points  of  shock-shock  interactions.  A 
wide  variation  of  Knudsen  numbers  may  exist  in  the  same  simulation  -  in  the  example  of  Fig.  3  the  Knudsen 
numbers  vary  from  0.026  on  the  toroid  to  0.88  on  the  tether.  This  range  indicates  a  variation  from  continuum  to 
near  free  molecular  flow  and  requires  analyses  appropriate  for  these  domains  in  the  same  simulation.  [21]  The 
inflation  infrastructure  provides  opportunity  to  engage  film  cooling  and  raise  acceptable  aerothermal  loads  with 
innovative  material  concepts. [22] 

2.1.3  Energetics 

Harvesting  and  storing  energy  imparted  to  the  gas  during  atmospheric  entry  (regenerative  aerobraking)  offers 
significant  potential  to  do  useful  work  after  landing.  Macheret  et.  al.  [23]  describe  a  concept  for  a  regenerative 
aerobraking  system  at  Mars  using  a  magnetohydrodynamic  (MHD)  system.  Though  not  usually  considered  an 
integral  part  of  planetary  vehicle  aerothermodynamics  this  potential  application  illustrates  that  MHD  simulation 
should  not  be  ignored  in  establishing  a  complete  entry  vehicle  design  competency. 

2.2  Exploration  of  Outer  Planets 

A  NASA  systems  study  for  aerocapture  at  Neptune  reveals  significant  challenges  to  the  thermal  protection 
system  design.  [24,  25]  The  study  concluded  that  aerocapture  at  Neptune  with  inertial  entry  velocity  of  29  km/s 
could  be  accomplished  with  a  rigid  aeroshell  configured  as  a  flattened  ellipsled  with  a  lift-to-drag  (L/D)  ratio 
of  0.80  and  a  ballistic  coefficient  j3  ~  895  kg/m2.  Laub  et.  al.  [25]  note  that  the  large  convective  and  radiative 
heating  rates  on  the  nose  (10-15  kW/cm2)  and  on  acreage  windside  (3-6  kW/cm2)  using  existing  materials 
require  thicknesses  that  currently  cannot  be  manufactured  as  a  high  quality  product.  Neither  the  effects  of 
coupled  ablation  and  radiation  nor  the  effects  of  significant  shape  change  on  aerodynamics  were  considered  in 
this  high  level  study  but  would  be  required  as  designs  mature  -  a  goal  of  CA  development  discussed  below. 

3.0  ALGORITHMS 

In  the  previous  chapter,  Problem  Definition,  we  discussed  simulation  capabilities  required  for  defining  the 
aerothermal  environments  of  planetary  entry  vehicles  and  opportunities  for  exploring  innovative  concepts  for 
EDL  systems.  In  this  chapter  we  highlight  two  algorithmic  topics  that  must  be  addressed  to  achieve  these 
simulation  capabilities.  We  note  these  topics  are  of  special  concern,  if  not  unique,  to  hypersonic  flow  simulation. 
The  first  addresses  the  need  to  obtain  accurate  heating  simulation  on  grid  systems  that  can  easily  adapted  to 
complex  flow  topologies.  The  second  addresses  the  infrastructure  needs  to  support  fully  coupled  ablation, 
radiation,  and  shape  change  in  a  CA  simulation.  Governing  equations  and  physical  models  for  atmospheric 
entry  have  been  presented  previously  [26]  and  will  not  be  repeated  here. 

3.1  Heating  on  Unstructured  Grids  in  Hypersonic  Simulations  —  Multi-Dimensional  Recon¬ 
struction 

The  simulation  of  hypersonic  flows  with  fully  unstructured  (tetrahedral)  grids  has  severe  problems  with  respect 
to  the  prediction  of  stagnation  region  heating.  [27,  28]  The  problems  arise  for  two  reasons. [29] 

First,  good  shock  capturing  in  the  hypersonic  regime  requires  alignment  of  the  grid  with  the  captured 
shock.  Alignment  here  means  that  all  surfaces  of  a  control  volume  in  contact  with  the  shock  are  either  parallel 
or  orthogonal  to  the  discontinuity.  Any  skewness  present  in  a  conventional  upwind  scheme  (in  which  first-order 
reconstruction  is  a  function  of  only  two  nodes  on  opposite  sides  of  a  shared  face)  results  in  a  smearing  of  the 
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flow,  which  manifests  itself  as  a  non-physical  distribution  of  entropy  from  streamline  to  streamline  crossing 
the  shock.  In  the  stagnation  region,  there  is  very  little  dissipation  of  these  entropy  gradients  as  they  pass  from 
the  shock  to  the  boundary-layer  edge,  and  the  surface  heating  is  very  sensitive  to  these  non-physical  gradients. 
Tetrahedral  cell  topologies  by  their  very  nature  will  always  have  at  least  one  surface  that  is  skewed  to  the 
captured  shock. 

Second,  tetrahedral  grids  are  not  well  suited  to  the  preservation  of  symmetry  and  the  biased  tetrahedral 
grids  used  in  the  simulations  of  heating  herein  are  particularly  ill-suited  for  this  purpose.  Asymmetric  elements 
and  faces  skewed  to  the  symmetry  surface  tend  to  induce  non-physical  cross  flows.  Again,  this  problem  is 
particularly  evident  in  the  stagnation  region  where  the  physical  flow  velocities  themselves  are  small  and  the 
non-physical,  grid-induced  velocities  are  of  similar  magnitude.  Surface  heating  rates  are  a  sensitive  indicator  of 
these  problems.  The  combination  of  random  jaggedness  in  the  shock  capturing  process  as  well  as  non-physical 
cross  flow  velocities  cross-couple  to  produce  unacceptable  surface  heating  predictions. 

Problems  with  heating  can  be  overcome  if  one  uses  a  semi- structured  grid  (e.g.,  prisms)  across  the  boundary 
layer  and  adapts  the  grid  to  the  shock.  Nompelis  et.  al.[30]  show  excellent  heating  results  on  families  of  grids 
where  a  prismatic  grid  was  adapted  to  the  shock  and  acceptable  heating  results  when  an  unbiased  (random  face 
orientation),  tetrahedral  grid  was  adapted  to  the  shock.  Prismatic  elements  are  relatively  easy  to  generate  on 
blunt  body  geometries  and  their  superior  performance  with  respect  to  aeroheating  predictions  have  led  to  their 
use  as  standard  practice  in  unstructured  simulations.  There  is  no  better  alternative  with  todays  algorithms  than 
to  use  prismatic  elements  to  capture  the  boundary  layer  and  bow  shock. 

If  one  is  doing  a  hypersonic  simulation  without  any  free  shear  layers  or  internal  shocks  the  specialized 
application  of  prismatic  elements  at  the  body  and  bow  shock  is  perfectly  acceptable.  If  such  flow  structures 
exist  —  as  they  do  in  almost  all  interesting  problems  noted  in  the  previous  chapter  —  then  the  accuracy  of 
any  algorithm  must  be  questioned  when  it  ignores  special  topological  grid  requirements  where  the  features  are 
harder  to  resolve.  Consequently,  overcoming  these  accuracy  issues  on  tetrahedral  grids  is  considered  one  of  the 
highest  priority  tasks  for  achieving  simulation  for  innovation. 

A  three-dimensional  reconstruction  algorithm  to  address  these  challenges  was  originally  described  in  ref¬ 
erence  [29].  Some  details  of  the  algorithm  are  included  in  the  Appendix  for  convenience  of  the  reader.  The 
key  element  of  the  algorithm  is  that  inviscid  flux  reconstruction  is  element-based  rather  than  edge-based  so  that 
even  the  first-order  reconstruction  is  based  on  4  nodes  surrounding  the  element  rather  than  2  nodes  defining  an 
edge.  We  note  that  this  cylinder  challenge  problem  and  a  related  sphere  challenge  problem  (not  included  here) 
have  been  presented  in  [31]. 


3.1.1  Challenge  Problem  -  Biased,  Tetrahedral  Grid  on  Cylinder  with  Spanwise  Resolution 

A  challenge  problem  is  identified  in  which  a  highly  biased,  tetrahedral  grid  is  applied  to  the  simulation  of  a 
hypersonic  flow  over  a  3D  cylinder  at  5km/s,  0.001kg/m3  and  T ^  200K  in  air  (M^  « 

17). [28,  32]  (In  fact,  other  hypersonic  free  stream  conditions  are  acceptable.  The  challenge  here  is  to  recover 
the  proper  spanwise  uniformity  within  1%.)  A  structured-grid  solution  generated  with  LAURA[33]  is  used 
as  both  a  benchmark  and  to  generate  initial  grids  for  use  in  FUN3D.[28]  A  successful  simulation  requires  a 
constant  spanwise  simulation  of  heating,  pressure,  and  shear  and  symmetric  simulation  of  these  quantities  on 
the  right  and  left  sides  of  the  cylinder.  The  grid  bias  accentuates  any  issues  with  the  algorithm.  To  date  there 
have  been  no  successful  simulations  of  this  challenge  problem  on  the  given  grid  using  conventional  upwind 
formulations.  (A  recent  thesis  on  a  related  unstructured  grid  without  bias  produced  excellent  heating  symmetry 
using  a  high-order  Discontinuous  Galerkin  finite  element  method  with  a  PDE  based  artificial  viscosity.  [34]) 
The  structured  grid  from  LAURA,  adapted  to  the  shock  and  boundary  layer,  is  converted  from  hexahedral 
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elements  to  tetrahedral  elements  by  adding  diagonal  edges  consistently  from  minimum  index  to  maximum  index 
corners.  A  comparison  of  the  grids  in  a  plane  of  nodes  perpendicular  to  the  cylinder  axis  is  presented  in  Fig.  4. 
The  structured  grid  has  65  nodes  from  the  body 
to  the  inflow  boundary  ahead  of  the  bow  shock 
and  61  nodes  from  the  left  to  right  outflow  bound¬ 
aries.  Node  placement  is  identical  between  the 
two  grids.  The  placement  of  additional  edges  in 
the  unstructured  grid  is  the  only  difference.  The 
strong  biasing  of  diagonals  in  the  grid  is  an  inten¬ 
tional  characteristic  to  expose  algorithm  deficien¬ 
cies  that  may  otherwise  be  averaged  out. 

A  key  element  of  this  test  problem  is  the  addi¬ 
tion  of  ten  spanwise  cells,  shown  in  Fig.  5,  across 
the  cylinder,  providing  additional  degrees  of  free¬ 
dom  in  the  simulation  to  allow  asymmetries  to  de¬ 
velop.  The  spanwise  degrees  of  freedom  enable 
non-physical  cross-flow  velocities  to  develop  and 
irregular  shock  capture  in  the  spanwise  direction 
that  combine  to  corrupt  the  predicted  heating  dis¬ 
tribution.  Earlier  tests  had  shown  that  single  span- 
wise  cell  simulations  were  not  severely  compromised. 


FUN3D 


LAURA 


Figure  4:  Structured  grid  (LAURA,  bottom)  and  biased,  unstruc¬ 
tured  grid  (FUN3D,  top)  in  plane  orthogonal  to  cylinder  surface. 


(a)  Structured,  quadrilateral  surface  elements 


(b)  Unstructured,  triangular  surface  elements 


Figure  5:  Surface  mesh  on  cylinder  with  ten  rows  of  spanwise  cells. 


A  large  spanwise  variation  in  heating  and  shear  for  this  simulation  using  conventional  reconstruction  is 
evident  in  Fig.  6(a)  with  ±30%  about  the  mean  —  an  unacceptable  result  for  aerothermodynamic  analyses. 
Recall  that  the  symbols  at  every  x  location  should  overplot  and  the  peak  value  of  shear  on  the  right  should  equal 
the  peak  value  of  shear  on  the  left.  The  LAURA  result  in  this  case  is  52  W/cm2  which  is  in  good  agreement 
with  the  spanwise  mean  of  the  heating  at  the  stagnation  point.  The  surface  heating  contours  in  Fig.  6(b)  indicate 
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that  a  slight  drift  in  the  velocity  —  probably  induced  by  grid  bias  -  produces  a  higher  heating  toward  the  front 
(//  =  0)  plane. 


(a)  Heating  (blue  squares)  and  shear  (red  squares)  distribution 
for  each  of  ten  rows  of  span  wise  cells. 


(b)  Surface  heating  contours. 


Figure  6:  Simulation  results  using  conventional  reconstruction  on  the  tetrahedral  grid  with  ten  spanwise  cells. 


(a)  Values  of  heating  -  blue  circles,  shear  -  red  triangles,  and 
pressure  -  black  squares  at  nodes  on  surface.  Figure  shows  results 
for  all  1 1  surface  nodes  at  a  given  x  location. 


q,  W/cm2:  10  15  20  25  30  35  40  45  50  55  60 


(b)  Heating  contours  on  surface  showing  significant  improve¬ 
ment  in  spanwise  symmetry.  Some  edge  effects  at  spanwise 
boundaries  still  persist. 


Figure  7:  Simulation  results  using  Option  2  reconstruction  on  the  tetrahedral  grid  with  ten  spanwise  cells. 


Results  for  the  challenge  problem  conditions  are  repeated  in  Fig.  7  using  multi-dimensional  reconstruction. 
The  new  results  (using  Option  2  as  defined  in  the  Appendix)  are  presented  as  symbols  and  a  benchmark, 
grid-converged,  structured  grid  solution  (LAURA)  for  the  same  case  are  presented  as  solid  lines.  The  heating 
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contours  show  excellent  spanwise  uniformity  in  Fig.  7(b)  in  the  interior.  Some  end  effects  are  evident  at  the  top 
and  bottom  boundaries.  The  heating  and  shear  distributions  are  in  good  agreement  with  the  benchmarks. 

The  element-based  reconstruction  algorithm  requires  three  reconstructions  per  element  whereas  the  base¬ 
line  ID  algorithm  requires  one  reconstruction  per  edge.  The  cylinder  challenge  problem  with  10  spanwise 
cells  required  14.8  seconds  per  relaxation  step  as  compared  to  6.8  seconds  per  step  for  the  baseline  ID  algo¬ 
rithm.  Some  improvement  is  expected  when  the  inviscid  and  viscous  terms  are  computed  in  a  single  loop  over 
elements.  This  refactoring  will  eliminate  some  redundancies  occurring  in  both  the  inviscid  and  viscous  formu¬ 
lations.  Additional  simulations  using  multi-dimensional  reconstruction  will  be  presented  in  the  next  chapter  on 
more  challenging  problems. 

3.2  Infrastructure  —  Free  Energy  Minimization  -  (FEM) 

The  Gibbs  Free  Energy  minimization  (FEM)  option  is  implemented  with  a  goal  of  making  simple  adjustments 
to  the  baseline  thermochemical  nonequilibrium  flow  solver  to  enable  a  thermochemical  equilibrium  solver  for 
generic  mixtures  of  gases  with  spatially  varying  elemental  mass  fractions.  The  capability  is  required  to  enable 
interfaces  with  ablation  routines  based  on  equilibrium  gas  chemistry.  The  capability  also  enables  simulations 
that  are  so  deep  into  the  equilibrium  regime  that  the  magnitude  of  species  chemical  source  terms  overwhelm 
the  magnitude  of  convection  and  diffusion.  (Roundoff  error  in  the  evaluation  of  the  source  terms  is  significant 
relative  to  the  magnitude  of  convective  and  diffusive  terms.)  Consequently,  the  sum  of  the  species  continuity 
equations  does  not  telescope  to  a  meaningful  mixture  continuity  equation  because  of  roundoff  error.  The  mod¬ 
ifications  are  originally  described  in  [35].  In  essence,  the  solution  of  Ns  species  continuity  equations  with  the 
solution  of  Ne  elemental  continuity  equations  and  ( Ns  —  Ne )  algebraic,  partial  equilibrium  equations. 


3.2.1  Elemental  Continuity  Equations 

The  transformation  from  species  continuity  to  elemental  continuity  involves  the  identification  of  the  fraction 
of  element  i  in  each  species  j.  The  transformation  matrix  F  involves  Ne  rows,  one  for  each  element,  and 
Ns  columns,  one  for  each  species.  The  fraction  is  computed  by  multiplying  the  elemental  weight,  by  the 
number  of  atoms  of  that  element  in  the  species,  and  dividing  by  the  species  weight  Mj. 


F  — 

rh3  ~ 


(1) 


Consider  the  one-dimensional  species  conservation  equations  set. 

dq  df 

- ( - —  w 

dt  dx 


(2) 


where  f  includes  convection  and  diffusion.  If  one  multiplies  both  sides  of  Eq.  2  by  F  one  obtains  the  elemental 
conservation  equations 


dq  df 

dt  +  dx 

dq  df 

dt  +  dx 


=  Fw 


0 


(3) 


(4) 


Note  that  Fw  =  0.  Chemical  reactions  cannot  create  or  destroy  an  element  in  the  control  volume.  Conse¬ 
quently,  there  is  no  need  to  add  the  source  term  to  the  species  continuity  equations  prior  to  assembling  the 
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elemental  continuity  equations.  Besides  saving  operation  counts  it  is  important  to  apply  the  transformation 
matrix  before  adding  the  source  term  because  the  critical  components  of  continuity  through  convection  and 
diffusion  may  be  of  the  order  of  roundoff  error  in  the  assembly  of  the  source  term  w.  The  assembly  of  the 
Jacobian  of  the  elemental  conservation  laws  (Eq.  12)  employs  the  same  matrix  transform  and  retains  the  use  of 
species  densities  (rather  than  elemental  densities)  as  the  fundamental  dependent  variables.  Consequently,  the 
remaining  non-base  species  must  be  solved  using  partial  equilibrium  relations  in  order  to  close  the  equation  set 
for  species  densities. 


3.2.2  Partial  Equilibrium  Relations 

The  derivation  of  partial  equilibrium  relations  through  the  minimization  of  free  energy  is  available  in  several 
texts.  [36]  Consider  the  partial  equilibrium  relation  for  non-base  species  j  which  may  be  expressed  as 

reactants  products 

T,  <=>  y  Pn,jNn  (5) 

m  n 


where  one  of  the  products  Nn  corresponds  to  non-base  species  j,  the  reactants  correspond  to  constituent  base 
species,  and  one  additional  species  may  be  required  to  enforce  a  chemical  balance.  The  partial  equilibrium 
relation  involving  non-base  species  j  may  then  be  derived  as  follows. 


AF» 


y!  Pn,j  ( 


Hn  -  TSn 
RT 


A  Ff  ^  .  RT  ^  .  RT 
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^  ^  fim,j  Pm  ^  ^  @n,j  Pnd ) 

m  n 

^  ^  /^m,  j  In  Pm  ^  ^  /5n,j  In  Pnd) 
m  n 

species 

^2  GjylnPk  (9) 

k 


where  R  =  pooJ?/101, 325  to  include  conversion  from  pressure  in  atmospheres  to  the  metric  system  and  to 
non-dimensionalize  density  by  The  non-base  species  continuity  equations  are  thus  replaced  by  a  linear 
equation  in  In  pp  with  coefficients  defined 


{0  if  species  k  is  not  part  of  reaction  j 

—f3nd  if  species  k  —  n  is  a  product  in  reaction  j  (10) 

/3mj  if  species  k  —  m  is  a  reactant  in  reaction  j 

The  temperature  dependent  equilibrium  constant  Kj  (T)  is  curve  fit  using 

Kj  (T)  =  A:J  ^  +  Bj  +  Cj  In  Z  +  D  j  Z  +  EjZ2  (1 1) 

where  Z  =  10, 000/T  and  the  constants  Aj  -  Ej  are  defined  by  fitting  to  the  left-hand- side  of  Eq.  8  (defining 
Kj(T )  in  the  left  hand  side  of  Eq.  9)  at  T  —  1000,  2000,  4000,  8000,  16,  000.  The  enthalpy  Hn  and  entropy 
Sn  for  species  n  is  obtained  from  curve  fits.  [37] 
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The  coupled  solution  of  Eq.  9  with  the  other  conservation  laws  can  be  a  challenge,  particularly  when  it 
is  early  in  the  simulation  and  the  solution  is  far  from  a  converged  state.  It  has  been  found  that  ignoring  the 
Jacobian  of  Eq.  9  with  respect  to  temperature  and  by  multiplying  the  Jacobian  of  Eq.  9  with  respect  to  species 
densities  by  a  safety  factor  of  15  has  enabled  solution  of  the  FEM  coupled  equation  set  for  all  tests  conducted 
for  this  paper.  The  factor  15  was  determined  through  numerical  experiment. 


3.3  Infrastructure  —  Ablation  Boundary  Conditions 

The  ablation  boundary  condition  derived  here  assumes  there  will  be  coupling  with  an  ablation  module  such  as 
FIAT[38]  in  conjunction  with  an  equilibrium  chemistry  module.  It  is  assumed  that  the  ablation  mass  flow  rate 
m  =  rhpyroiysis  +  rn char  ^  the  elemental  composition  of  ablation  gases  and  the  temperature  of  the  ablating 
surface  are  supplied  by  the  ablation  module  or  can  be  user  specified.  The  ablation  module  itself  requires  the 
aerothermodynamic  environments  (pressure,  radiative  and  convective  heating,  and  possibly  shear)  over  time  to 
compute  material  response.  Consequently,  a  boundary  condition  must  be  applied  to  the  flow  solver  that  closes 
the  equation  sets  for  either  nonequilibrium  external  flow  using  species  continuity  equations  or  equilibrium 
external  flow  using  elemental  continuity  equations.  As  noted  in  the  previous  section  the  species  densities  are 
retained  as  the  fundamental  dependent  variable  data  set  regardless  of  the  assumption  of  equilibrium  state  of  the 
external  flow. 

It  is  assumed  that  the  gas  chemistry  at  and  within  the  ablating  heatshield  is  in  an  equilibrium  state.  A  steady, 
one-dimensional,  elemental  continuity  equation  (Eq.  12  )  describes  the  flow  of  an  equilibrium  gas  across  the 
boundary  that  accounts  for  convection  and  diffusion  of  element  i. 


dpjv  _  dpjVj 
dy  dy 


(12) 


Because  we  wish  to  work  directly  with  species  densities  the  transformation  matrix  introduced  in  the  previous 
section  is  applied  in  Eq.  13  below. 


dpjv 

J  dy 


=  F, 


h3 


dpjVj 

dy 


(13) 


where  summation  over  species  j  is  implied. 

Integrate  Eq.  13  over  y  from  some  virtual  point  sufficiently  deep  in  the  heatshield  where  it  is  assumed  that 
the  net  elemental  mass  flux  is  given  by  c^/ra.  This  assumption  ignores  the  possibility  that  ablating  mass  flux 
sources  are  distributed  through  the  depth  of  the  char.  Short  of  simulating  flow  through  the  char  it  is  not  clear 
that  this  modeling  assumption  can  be  improved.  The  resulting  relation  (Eq.  14)  provides  a  balance  for  each 
element  in  which  the  convected  elemental  mass  flux  out  of  the  boundary  minus  the  diffused  elemental  mass 
flux  into  the  boundary  must  equal  the  ablated  elemental  mass  flux  specified  for  the  material. 


FijpjV  -  Citabim  =  FijpjVj  (14) 

F i,jPjV  =  FijpDj-^—  (15) 

Note  that  an  effective  binary  diffusion  coefficient  is  applied  to  model  diffusive  flux.  A  mass  corrected  diffu¬ 
sive  flux[39]  is  also  applied  to  guarantee  that  the  sum  of  diffusive  flux  over  all  species  (or  over  all  elements) 
goes  to  zero.  Simulations  show  some  elemental  separation  near  the  wall  for  the  no  blowing  case.  Elemental 
separation  in  the  no-blowing  case  is  less  when  modeling  multi-component  diffusion  using  the  Stefan  Maxwell 
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equations.  [39]  A  verification  check  shows  that  no  elemental  separation  is  computed  in  the  no-blowing  case 
when  a  constant  Schmidt  number  is  applied. 

The  discretization  of  Eq.  15  is  shown  in  Eq.  16.  It  is  necessary  to  solve  for  all  dependent  variables  in  the 
ghost  cell  0. 


h3 


(pjv)  1  +  (pjV)  0 


C itabl —  Fi,j 


(pDj)i  +  (pDj)  o 


Xj,  1  ~  Xj, 0 

Ay 


(16) 


The  solution  in  the  ghost  cell  requires  simultaneous  relaxation  of  the  partial  equilibrium  relations  given  by 
Eq.  9.  The  Jacobian  formed  using  In  p3  as  the  dependent  variable  provides  a  more  robust  convergence  of  this 
boundary  condition.  A  relaxation  factor  is  applied  to  each  row  of  the  Jacobian  according  to  the  magnitude  of  the 
current  residual  for  that  row  equal  to  n/(d  In [pj]2+n)  where  n  is  the  local  iteration  count  on  relaxation  steps  for 
this  boundary  condition.  Note  that  in  the  limit  of  zero  blowing  this  boundary  condition  specifies  an  equilibrium 
catalytic  wall  -  that  is,  the  species  are  in  chemical  equilibrium  as  determined  by  the  wall  temperature  and 
pressure  and  elemental  mass  fraction  gradient  equal  to  zero. 


1 

2 


[(Hi  +  (Ho] 


=  m 


(17) 


Pi  +  pivj=  PO  +  poVo  (18) 

The  equation  set  is  closed  with  a  discretization  of  net  mass  flux  (Eq.  17)  and  normal  momentum  balance 
(Eq.  18).  Eq.  17  is  used  to  eliminate  vq  from  all  equations.  Then  Eq.  18  replaces  one  of  the  elemental  continuity 
equations  so  that  the  ghost  cell  pressure  po  is  consistent  with  the  specified  temperature  Tq  =  2 Tw  —  T\  and  the 
computed  densities  p3$. 

At  present,  ablation  rates  and  ablation  gas  elemental  mass  fractions  are  explicitly  set  by  the  user.  The 
boundary  condition  is  sufficiently  robust  to  enable  a  coupled  simulation  assuming  equilibrium  gas  chemistry  at 
Voo  =15  km/s  and  a  dimensionless  blowing  rate  equal  to  0.2. 

On  the  first  step  of  initialization,  all  cells  (including  ghost  cells)  are  set  to  inflow  conditions.  On  the  second 
step  of  initialization  ghost  cell  velocities  are  reset  to  zero  and  wall  temperatures  are  set  to  a  user  specified 
value.  Simulations  involving  ablation  and/or  radiation  are  started  from  a  converged  (or  nearly  converged)  non¬ 
radiating  and  non-ablating  solution.  The  line-implicit  formulation  for  relaxation  of  interior  cells  requires  a 
Jacobian  of  the  boundary  condition.  It  is  assumed  for  this  purpose  that  species  mass  fractions  and  temperature 
in  the  ghost  cells  are  fixed.  (The  point-implicit  formulation  for  relaxation  of  interior  points,  by  definition, 
requires  no  Jacobian  information  of  neighbor  cells.) 

During  interior  cell  updates  the  surface  convective  heating,  shear,  and  diffusion  flux  are  calculated  from 
current  values  in  the  ghost  cells.  Ghost  cell  values  in  turn  are  updated  at  a  user  defined  frequency,  Nexchange. 
(Ghost  cells  are  also  used  to  hold  information  from  neighbor  blocks  and  the  variable  name  Nexchange  derives 
from  the  frequency  in  which  neighbor  block  information  is  exchanged.)  Ghost  cell  values  of  species  density 
are  fully  defined  by  interior  point  values  -  there  is  no  time  dependent  update  or  damping  using  any  fraction  of 
a  previous  value.  Damping  is  applied  to  the  update  of  a  radiative  equilibrium  surface  temperature  in  the  case 
of  no  blowing.  The  radiative  equilibrium  wall  temperature  boundary  condition  is  then  defined  with  Twnew  = 
e[q/ (e<r)] 1/4  +  (1  —  e)Twold .  In  this  case,  the  ghost  cell  temperature  is  computed  as  Tq  =  2 Tw  —  T\. 

N exchange  is  set  equal  to  2  in  all  cases  presented  herein.  The  temperature  damping  factor  e  is  varied  between 
0.001  and  0.01.  In  more  recent  applications  (some  involving  a  surface  energy  balance  not  covered  here)  the 
combination  Nexchange  —  20  and  s  —  0.1  have  provided  more  robust  convergence. 
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3.4  Infrastructure  —  Radiation  Coupling 

Radiative  energy  transport  is  computed  using  the  tangent  slab  approximation  with  the  HARA  code.  [40,  41]. 
The  flow  solver  provides  species  densities,  temperatures,  and  distances  along  a  ray  orthogonal  to  the  surface 
as  inputs  to  HARA.  HARA  returns  divergence  of  the  radiative  flux  and  radiative  flux  to  the  wall  for  use  in 
LAURA.  Calls  are  made  to  HARA  on  the  order  of  hundreds  of  relaxation  steps,  Nrad,  holding  the  radiative 
quantities  constant  over  this  period  (loosely  coupled).  Options  are  provided  to  interpolate  radiative  quantities 
between  selected  rays  to  reduce  overall  time  in  the  radiation  solver.  The  HARA  radiation  model  treats  atomic 
lines  and  continuum  of  N,  O,  C,  and  H.  Molecular  bands  are  treated  for  N2,  N^,  O2,  NO,  C2,  C3,  CN,  CO,  C2H 
with  the  smeared-rotational  band  (SRB)  method  for  molecular  band  radiation.  Non-Boltzmann  electronic  state 
models  are  applied  for  atoms  and  molecules.  Typical  computational  time  for  a  single  line-of-sight  is  20  to  50 
seconds  -  short  enough  to  enable  fully  coupled  simulations  with  codes  like  LAURA. 

A  parameter  frac[  may  be  specified  to  introduce  the  divergence  of  the  radiative  flux  more  slowly  into  the 
flow  solver. 


V  •  <ad  =  fradV  ■  qrad(HARA)  +  (1  -  /rad)V  •  q?ad  (19) 

If  the  radiation  introduces  strong  cooling  into  the  shock  layer  the  shock  standoff  distance  decreases  signif¬ 
icantly.  As  the  captured  shock  moves  closer  to  the  body,  computational  cells  previously  in  the  hot  shock  layer 
are  overtaken  by  the  pre-shock  domain.  However,  these  cells  retain  the  cooling  term  V  •  q rad,i,j,k  from  the 
previous  call  to  the  radiation  solver.  The  LAURA  code  enforces  a  floor  on  temperature  for  every  point  in  the 
flowfield  that  allows  the  simulation  to  survive  such  transients  under  most  circumstances.  However,  cases  have 
been  observed  where  introduction  of  the  full  update  to  radiation  has  caused  stability  problems.  Best  practice  so 
far  has  been  to  set  frad  to  0.4  on  the  first  call  from  an  uncoupled  solution  and  then  to  increase  its  value  to  1  on 
subsequent  returns  from  the  radiation  solver. 

4.0  SIMULATIONS  FOR  VALIDATION 
4.1  Validation  Challenges 

Let  us  focus  on  the  unique  challenges  of  code  validation  for  simulation  of  hypersonic  flow  for  planetary  entry. 
The  key  issue  in  the  hypersonic  flight  domain  is  to  adequately  model  the  mechanisms  through  which  a  massive 
reservoir  of  kinetic  energy  cascades  across  diverse  chemical  and  internal  modes.  More  precisely,  the  kinetic 
energy  of  an  entry  vehicle  cascades  into  chemical,  translational,  and  internal  modes  of  the  quiescent  atmospheric 
gas  through  action  of  the  vehicle’s  bow  shock.  The  transfer  of  energy  through  space  by  laminar  and  turbulent 
transport,  convection  and  radiation  are  related  critical  elements  in  a  standard  simulation.  Add  to  this  mix  the 
physics  of  gas-surface  interactions  (catalysis),  material  response  (ablation),  jet  interactions  (retro-propulsion), 
and/or  elasto-dynamic  boundary  conditions  (inflatable  /  deployable  structure)  and  the  challenges  to  planetary 
entry  validation  are  daunting. 

The  path  forward  must  combine  measured  data  from  ground-based  tests  and  flight  experiments  as  well  as 
computed  data  from  first-principles  physics  simulations  on  canonical  problems  (computational  chemistry [42], 
DNS[43]).  Each  has  pros  and  cons.  We  confine  our  attention  here  to  measured  data. 

Ground-based  tests  provide  the  best  opportunity  to  measure  boundary  conditions  and  probe  the  flow  field 
with  sufficient  detail  to  validate  physical  models.  However,  the  facilities  are  not  generally  quiet  and  often  the 
flows  are  not  perfectly  uniform  nor  steady.  As  enthalpies  are  increased  flow  quality  and  test  times  become 
even  more  problematic.  Recent  efforts  by  Candler  et  al.  [44]  underscore  the  proposition  that  future  code 
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validation  tests  will  require  more  complete  simulation  of  the  facility  flow  in  addition  to  simulation  around 
models.  Simulation  of  facility  flow  with  requisite  validation  checks  should  become  a  standard  component 
of  CA  analysis  practice.  [45]  One  thus  acknowledges  that  ground-based  facilities  often  “bruise”  the  gas  (flow 
uniformity,  chemical  and  thermal  state)  but  in  fact  these  perturbations  provide  additional  opportunity  to  validate 
the  suite  of  physical  models  used  in  CA  simulations. 

As  noted  by  Wright  el.  al.  [46]  “Flight  testing,  while  certainly  capable  of  fully  exercising  the  TPS  design, 
is  prohibitively  expensive  in  most  cases.  In  general  flight  testing  should  be  reserved  for  final  model  and  system- 
level  validation,  once  we  have  good  physics-based  models  of  the  expected  environment  and  resulting  material 
performance.  Ideally,  the  purpose  of  a  dedicated  flight  experiment  should  not  be  to  learn  anything  new  about  the 
flight  environment  encountered,  but  rather  to  validate  that  the  models  employed  are  adequate.”  However,  valu¬ 
able  engineering  information  can  be  obtained  at  relatively  lower  cost  by  including  engineering  instrumentation 
on  science  missions.  While  there  is  little  benefit  to  the  mission  that  flies  it,  the  overall  exploration  program  ben¬ 
efits  from  the  additional  knowledge  gained  about  the  flight  environment.  This  new  knowledge  may  then  benefit 
future  missions  through  reduction  of  TPS  mass  requirements,  for  example,  thus  enabling  increased  payload 
mass.  The  Mars  Science  Laboratory  Entry,  Descent,  and  Landing  Instrumentation  (MEDLI)  Project  is  a  prime 
example  of  piggybacking  engineering  instrumentation  on  a  science  mission  -  Mars  Science  Laboratory.  [47] 
The  flight  science  objectives  of  MEDLI  directly  address  the  largest  uncertainties  in  the  ability  to  design  and 
validate  a  robust  Mars  entry  system,  including  aerothermal,  aerodynamic  and  atmosphere  models,  and  thermal 
protection  system  (TPS)  design.  [47] 

The  remainder  of  this  chapter  provides  example  of  code  assessment  and  validation  using  the  structured  flow 
solver  LAURA  and  the  unstructured  flow  solver  FUN3D  in  problems  of  particular  interest  to  the  CA  community. 
Both  ground-based  and  flight  data  are  simulated.  As  CA  practice  matures  the  simulation  of  each  of  these  cases 
should  become  routinely  accessible  with  published  uncertainties  of  predicted  aerothermal  environments. 

4.2  Retro-Propulsion 

Simulations  are  executed  to  assess  the  ability  of  CFD  to  predict  the  flow  field  associated  with  a  jet  firing  from 
a  blunt  body  into  on  oncoming  supersonic  flow.  Of  particular  interest  is  the  ability  to  achieve  a  grid  converged 
solution  in  a  flow  topology  where  the  location  of  critical  domains  requiring  fine  resolution  are  not  known  apriori. 
The  sensitivity  of  flow  topology  and  aerothermal  loads  as  a  function  of  turbulence  modeling  are  also  quantified 
in  these  studies.  The  test  case  involves  a  Mach  2  oncoming  flow.  ( —  527m/s,  —  3.55  x  10-2kg/m3) 

[48].  Experimental  data  indicates  this  condition  results  in  a  steady  interaction  -  removing  one  significant  level 
of  complexity  from  the  simulation. 

The  simulations  in  Fig.  8  used  the  unstructured  grid  solver  FUN3D.  The  initial  grid  generation  process  was 
found  to  be  critical  in  this  application.  Grid  adaptation  driven  by  local  gradient  information  failed  to  produce 
grid  converged  solutions.  Rather,  the  grid  converged  on  intermediate  flow  topologies  and  then  locked  in  on 
these  states  -  failing  to  evolve  through  a  coarser  grid  that  bounds  the  high  gradient  regions  of  the  intermediate 
states.  Ultimately,  the  expert  user  was  required  to  throw  a  relatively  fine  initial  grid  across  the  entire  interaction 
domain  with  judicious  use  of  source  terms  in  the  grid  generation  procedure  before  a  grid  converged  solution 
with  qualitative  agreement  with  experimental  data  could  be  obtained.  (This  situation  is  repeated  in  the  next 
test  problem  -  the  sharp  double  cone.  A  fine  initial  grid  is  required  to  enable  the  separation  zone  to  progress 
upstream  as  far  as  it  needs.)  Approximately  4.3  million  nodes  (25  million  tetrahedra)  were  used  to  resolve  one 
quarter  of  the  flow  domain  surrounding  the  model  (zero  degrees  to  ninety  degrees  around  the  symmetry  axis). 
The  requirement  here  for  an  expert  in  the  loop  (in  this  case  an  expert  in  training!)  to  recognize  that  the  grid 
adaptation  strategy  was  converging  on  a  wrong  solution  is  a  show  stopper  to  simulation  for  innovation.  How 
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(a)  Turbulent  shear  layer  -  SST  model 

-3 


(b)  Turbulent  shear  layer  -  Spalart  -  Almaras  model 

o  1 


(c)  Laminar  flow 

Figure  8:  Mach  contours  in  simulation  of  retro-propulsion  in  a  Mach  2.5  oncoming  flow. 


can  one  explore  innovative  concepts  for  thruster  configurations  when  the  simulation  for  any  one  configuration 
requires  so  much  oversight?  The  application  of  adjoint-based  adaptation[49-51]  is  expected  to  overcome  this 
problem  of  convergence  to  the  wrong  answer  and  thus  is  a  high  priority  task  on  the  path  to  simulation  for 
innovation. 

With  the  requisite  initial  grid  generation  strategy  established  it  was  possible  to  study  effects  of  turbulence 
modeling  on  the  interaction.  These  scoping  results  using  the  SST  model  [52]  (Fig.  8(a)),  the  Spalart  Almaras 
model  [53]  (Fig.  8(b)),  and  laminar  flow  (Fig.  8(c))  show  significant  variation  of  flow  topology.  Indeed,  the 
variation  is  not  unlike  that  seen  with  the  different  initial  grids  discussed  above.  Results  using  the  SST  model 
are  in  closest  qualitative  agreement  with  the  experimental  data.  Significance  of  agreement  is  only  qualitative 
because  only  Schlieren  results  were  available  and  some  inflow  conditions  were  not  defined.  A  new  experimental 
validation  program  has  been  initiated  and  supporting  simulations  are  being  used  to  assess  possibility  of  flow 
blockage. 
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4.3  Sharp,  Double  Cone 


A  series  of  experiments  to  measure  pressure  and  heating  for  code 
validation  involving  hypersonic,  laminar,  separated  flows  over  dou¬ 
ble  cones  was  conducted  at  the  Calspan-University  at  Buffalo  Re¬ 
search  Center  (CUBRC)  in  the  Large  Energy  National  Shock  (LENS) 
tunnel.  [54]  The  configuration  (Fig.  9)  yields  a  flow  topology  with  a 
large  separation  zone  across  the  intersection  of  the  two  cones.  An 
overview  of  the  flow  topology  is  evident  in  the  temperature  and  pres¬ 
sure  contours  of  Fig.  10.  The  shock  off  the  forward  part  of  the  sepa¬ 
ration  intersects  the  bow  shock  over  the  second  cone.  A  transmitted 
shock  from  this  intersection  reflects  off  the  wall  to  a  contact  surface. 

The  extent  of  separation  and  the  location  of  pressure  maxima  down¬ 
stream  of  the  reattachment  point  are  sensitive  to  gas  chemistry.  These 
features  are  also  a  sensitive  function  of  grid  convergence.  This  sen¬ 
sitivity  to  physical  models  and  grid  convergence  makes  simulation  of 
this  flow  an  excellent  code  validation  test.  The  experiments  were  the 

focus  of  a  blind  code  validation  study.  [54]  A  single  case,  Run  28,  is  featured  here  with  inflow  conditions 
Voo  =  2.664  km/s,  =  6.55  10-4  kg/m3,  =  185.6  K,  and  wall  temperature  Tw  =  293.3  K.  These 
conditions  correspond  to  M ^  —  9.59  and  Re^  =  144, 010  /m  in  Nitrogen  (to  minimize  effects  of  chemical 
reactions).  The  conditions  are  the  same  as  originally  presented  to  the  study  participants.  A  subsequent  evalua¬ 
tion  of  test  conditions  was  spurred  by  a  consistent  disagreement  (20%)  between  experimental  data  for  heating 
on  the  front  cone  prior  to  separation  and  multiple,  grid  converged  simulations.  [55]  This  study  found  that:  (1) 
vibrational  freezing  of  flow  exiting  the  facility  nozzle;  (2)  vibrational  energy  accommodation  at  the  surface; 
and  (3)  small  levels  of  flow  non-uniformity  must  be  included  in  the  simulations  to  improve  the  heating  predic¬ 
tion.  These  corrections  to  boundary  conditions  are  not  included  in  these  results.  The  intent  here  is  to  show 
that  the  simulations  on  a  tetrahedral  grid  system  are  consistent  with  the  earlier,  grid  converged,  structured-grid 
simulations. 


Figure  9:  Dimensions  (cm)  of  sharp,  double 
cone. 


The  computed  pressure  coefficients  and  Stanton  number  distributions  from  three  simulations  are  compared 
to  experimental  data  in  Fig.  11.  All  three  simulations  use  the  same  set  of  nodes  and  Roe  /  STVD  based  invis- 
cid  flux  formulation.  The  LAURA  solution  was  generated  for  the  blind  validation  study  (solid  blue  line)[56] 
and  uses  a  cell-centered,  structured-grid  formulation.  The  FUN3D  simulation  on  hexahedral  grid  (small  black 
circles)  uses  conventional,  one-dimensional,  edge-based  inviscid  flux  formulation.  The  FUN3D  simulation 
on  tetrahedral  grid  (small  green  triangles)  uses  the  new  three-dimensional,  element-based  inviscid  flux  for¬ 
mulation.  The  “LAURA  Path”  descriptor  in  these  figures  indicated  that  FUN3D  uses  the  same  modules  for 
thermodynamics  and  transport  properties  as  the  original  LAURA  code.  The  3D  reconstruction  on  tets  indicates 
a  slight  increase  in  the  extent  of  separation.  The  FUN3D  hex  simulation  falls  almost  exactly  under  the  original 
LAURA  simulation.  The  primary  differences  between  the  two  hex  simulations  is  that  LAURA  is  cell-centered 
and  uses  a  directionally  dependent  entropy  fix  while  the  ID  reconstruction  used  for  hex  grids  in  node  cen¬ 
tered  FUN3D  uses  a  directionally  independent  entropy  fix.  The  peak  pressure  appears  beneath  the  transmitted 
shock  at  x/L  «  1.3.  All  three  simulations  begin  the  rise  to  peak  pressure  slightly  after  the  experimental  data 
(in  red  squares).  This  delay  is  indicative  of  the  slightly  larger  extent  of  separation  with  the  original  boundary 
conditions.  The  reflected  waves  in  all  three  simulations  closely  follow  each  other  downstream  of  the  peak  pres¬ 
sure  and  heating.  This  result  demonstrates  the  consistency  of  the  3D  reconstruction  on  tetrahedral  grids  with 
both  structured  and  unstructured  grid  simulations  using  hexahedral  grids  for  a  complex  flow  in  which  principal 
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(a)  Temperature  contours. 


(b)  Detailed  view  of  pressure  contours  and  streamlines  in  separa¬ 
tion  region. 


Figure  10:  Flowfield  over  sharp,  double  cone. 


x/L 


(a)  Pressure  coefficient.  (b)  Stanton  number.. 

Figure  11 :  Comparison  of  experimental  data  with  computation  for  the  sharp,  double  cone,  Run  28. 


directions  vary  significantly  across  the  domain. 

All  three  simulations  employed  local-time- stepping  (constant  CFL).  When  starting  from  a  uniform  flow 
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this  process  produces  a  fully  attached  flow  (no  separation)  with  shocks  set  close  to  the  surface  in  the  early 
stages  of  the  simulation.  As  the  relaxation  process  continues,  the  shocks  move  out  from  the  body  and  the 
primary  separation  bubble  forms  at  the  intersection  of  the  two  cones  and  slowly  grows.  In  both  the  LAURA 
and  FUN3D  hexahedral  grid  simulations  the  separation  bubble  continues  to  grow  until  it  reaches  its  converged 
state.  In  the  3D  reconstruction  on  tetrahedral  grids  the  growth  of  the  separation  bubble  followed  much  the  same 
history.  However,  when  the  leading  edge  of  the  separation  bubble  reached  x/L  «  0.7  the  bubble  shape  became 
unsteady,  giving  way  to  a  pulsing  that  could  not  be  quenched  with  any  variation  of  CFL  or  entropy  fix.  A  steady 
solution  could  only  be  recovered  by  engaging  a  time  accurate  simulation  with  local  sub-iterations  employed 
between  time  steps.  The  use  of  time  accuracy  had  been  required  by  LAURA  for  a  related  flow  condition  at  a 
higher  Reynolds  number  (Run  24)  of  the  original  study.  The  3D  reconstruction  appears  to  be  more  sensitive  to 
this  phenomena. 

Expert  oversight  requirements  for  this  problem  are  somewhat  counter-intuitive.  The  LAURA  code  has  two 
grid  adaptation  features  that  were  tested  on  the  initial  simulation.  The  first  feature  redistributes  points  across 
the  boundary  layer  and  shock  layer  according  to  a  user  specified  cell  Reynolds  number  criteria  and  percentage 
of  points  requested  in  the  boundary  layer.  This  first  feature  worked  flawlessly  on  this  simulation.  The  second 
feature  has  a  bow  shock  sensing  algorithm  that  adjusts  the  length  of  the  computational  line  of  points  (structured 
grid)  such  that  the  inflow  boundary  is  aligned  with  the  captured  bow  shock.  This  feature  actually  caused  more 
problems  with  the  evolution  of  the  simulation.  The  shock  transmitted  off  the  leading  edge  of  the  separation 
bubble  caused  significant  readjustment  of  the  outer  boundary.  Each  grid  update  induced  changes  in  the  separa¬ 
tion  point.  The  feedback  process  —  change  in  separation  point  to  change  in  outer  boundary  and  grid  position  to 
change  in  separation  point ...  —  sustained  large,  pulsing  motion  across  the  cone.  Ultimately,  the  best  solution 
was  to  apply  uniform  grid  refinement  across  the  entire  domain  until  a  grid  converged  solution  was  obtained  and 
forego  any  attempt  to  save  grid  resources  with  shock  alignment.  These  observations  suggest  this  problem  is  an 
excellent  test  bed  for  grid  adaptation  algorithms  ultimately  required  in  simulation  for  innovation.  The  uniform 
refinement  approach  required  a  grid  of  approximately  512  x  1000  cells  to  demonstrate  grid  convergence  for  the 
256  x  500  solution.  An  intriguing  challenge  is  to  execute  a  grid  converged  solution  with  robust  adaptation  using 
an  order  of  magnitude  fewer  points.  Here  again,  application  of  adjoint-based  adaptation  needs  to  be  proven  on 
such  a  challenge  problem. 


4.4  STS  -  2 

CFD  simulations  of  heating  over  the  shuttle  orbiter  (STS  -  2)  including  effects  of  finite  surface  catalysis  have 
been  reported  for  the  multi-block,  structured  grid  solver  LAURA.  [57]  Some  of  these  solutions  were  redone 
with  a  refined  grid  and  updated  physical  models  in  support  of  the  Columbia  accident  investigation  and  Re¬ 
turn  to  Flight  activities.  [58]  The  original  simulations  used  a  grid  with  approximately  1  million  cells  and  a 
7-species  air  model  with  an  early  model  for  surface  catalysis.  The  more  recent  simulations  on  a  refined  grid 
with  approximately  8  million  cells  used  a  more  current  catalysis  model[59]  for  reaction-cured  glass  (RCG)  and 
a  5-species  air  model  (NO+  and  e-  were  ignored  to  save  computational  time  with  focus  on  relative  effects 
due  to  damage).  Two  new  simulations  are  run  with  FUN3D  for  a  Mach  24  case[57]:  Voo  =  6.920  km/s, 
Poo  =  5.750  10-5  kg/m3,  T ^  =  202  K,  a  =  39.4  deg  with  wall  temperature  defined  assuming  radiative  equi¬ 
librium  with  a  surface  emissivity  of  0.89.  The  structured  grid  was  converted  to  FUN3D  format  and  a  simulation 
was  executed  using  edge-based  reconstruction  on  hexahedral  elements.  The  hexahedral  grid  was  converted  into 
a  pure  tetrahedral  system  with  identical  nodes  and  a  second  FUN3D  simulation  was  executed  —  this  time  using 
element-based  reconstruction  on  the  biased,  tetrahedral-grid  system. 
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Two  cuts  along  the  windside  of  the  shuttle  orbiter  are  illustrated  in 
Fig.  12  over  a  flooded  contour  plot  of  surface  heating.  Heating  levels 
from  the  original  LAURA  simulation  from  1993,  an  updated  LAURA[60] 
simulation,  and  the  FUN3D  simulation  on  a  hexahedral  grid  and  a  tetra¬ 
hedral  grid  are  compared  in  Fig.  13  along  these  cut  lines.  The  updated 
LAURA  solution  for  heating  is  in  excellent  agreement  with  the  FUN3D 
simulation  on  a  hexahedral  grid.  These  two  simulations  use  identical 
modules  for  thermodynamics,  transport  properties,  chemical  kinetics, 
and  surface  catalysis.  The  heating  distribution  with  these  new  simula¬ 
tions  are  generally  in  good  agreement  with  the  experimental  data.  These  simulations  did  not  include  a  deflected 
body  flap.  Though  not  shown  here,  comparisons  along  other  cuts  across  the  body  (windside,  leeside,  circum¬ 
ferential)  show  equivalent  trends.  The  older  LAURA  results  over-predict  both  the  flight  data  and  the  newer 
simulations  along  these  windside  cuts.  It  is  thought  that  the  updated  model  for  surface  catalysis  is  the  primary 
reason  for  this  difference  though  a  systematic  study  has  not  been  executed  to  check  this  supposition.  The  main 
point  of  this  comparison  is  to  demonstrate  that  the  node  based  FUN3D  provides  equivalent  results  to  the  cell 
based  LAURA  when  identical  grids  and  physical  models  are  employed. 


Figure  12:  Cut  lines  over  windside  of 
shuttle  orbiter  showing  surface  heating 
contours. 


STS-2  Mach  24.3  STS-2  Mach  24.3 


(a)  y=0  in. 


(b)  y=185.  in. 


Figure  13:  Comparison  of  heating  on  windside  shuttle  test  case  among  original  LAURA  (black  line),  updated  LAURA  version  5 
(green  line),  FUN3D  on  hex  grid(red  line),  FUN3D  on  tet  grid  with  3D  reconstruction  (blue  line),  and  experimental  data  (black  error 
bars). 


The  next  step  is  to  demonstrate  that  the  FUN3D  simulation  on  tetrahedral  grids  with  3D  reconstruction 
provides  equivalent  surface  heating  as  the  hexahedral  grid  simulation.  This  step  has  been  obstructed  by  a 
sensitivity  to  the  formulation  of  the  outflow  boundary  condition  off  the  trailing  edge  of  the  leeside  of  the  wing 
and  vertical  tail.  The  boundary  condition  has  been  revised  to  vacuum  type  boundary  (enforce  low  pressure  on 
exit)  if  flow  reversal  is  detected  through  the  outflow  boundary.  The  intent  is  to  nudge  a  reverse  flow  to  exit  the 
outflow  boundary.  The  simulation  runs  with  this  fix  though  some  irregularities  in  the  heating  near  the  outflow 
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boundary  on  the  leeside  wingtip  are  evident. 


Figure  14:  Comparison  of  surface  contours  on  windside  of  shuttle.  Upper  surface  to  symmetry  plane  was  computed  using  a 
hexahedral  grid  system.  The  lower  surface  to  the  symmetry  plane  was  computed  using  biased  tetrahedral  grid  system. 


Of  greater  concern  is  the  appearance  of  a  narrow 
band  of  low  heating  along  the  plane  of  symmetry.  In 
Fig.  14(a)  we  show  heating  contours  on  the  lower  sur¬ 
face  of  the  shuttle  orbiter.  The  upper  half  of  the  fig¬ 
ure  was  simulated  with  FUN3D  using  a  hexahedral 
grid  system.  The  lower  half  of  the  figure  was  sim¬ 
ulated  independently  with  a  biased,  tetrahedral  grid 
system  using  identical  nodes.  Ideally,  the  bottom  half 
should  reflect  the  top  half.  In  this  case,  heating  con¬ 
tours  near  the  symmetry  plane  on  the  lower,  tetrahe¬ 
dral  grid  side  abruptly  run  up  toward  the  nose  indicat¬ 
ing  a  cooler  region  near  the  symmetry  plane.  Further 
away  from  the  symmetry  plane  the  heating  contours 
on  top  and  bottom  appear  to  be  in  reasonably  agree¬ 
ment  as  confirmed  in  viewing  the  blue  line  plots  in 
Fig.  13(b).  Note  that  pressure  contours  in  Fig.  14(b) 
show  reasonable  reflective  properties  across  the  sym¬ 
metry  plane  -  indicating  again  that  heating  and  not 
pressure  should  be  used  as  a  metric  of  solution  qual¬ 
ity  in  hypersonic  simulations.  Heating  distibutions 
on  the  leeside  (Fig.  15),  while  not  perfectly  reflec¬ 
tive,  do  not  reveal  any  inherent  problems  in  execution 


Figure  15:  Comparison  of  surface  heating  contours  on  leeside  of 
shuttle.  Lower  surface  to  symmetry  plane  was  computed  using 
a  hexahedral  grid  system.  The  upper  surface  to  the  symmetry 
plane  was  computed  using  biased  tetrahedral  grid  system. 
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using  tetrahedral  grids.  Note  the  heating  irregularities  at  the  leeside  wingtip  reported  earlier. 

The  cause  of  the  cool  streak  down  the  windside  centerline  in  the  tetrahedral  simulation  is  unknown.  A 
suspected  cause  -  an  asymmetry  in  the  way  the  eigenvalue  limiting  is  engaged  as  a  function  of  principal  direction 
-  has  been  investigated  by  making  radical  changes  in  the  limiting  process.  Radical  changes  have  also  been  made 
in  the  definition  of  principal  direction.  In  all  cases,  the  quality  of  the  surface  contour  pattern  is  not  significantly 
altered.  Future  tests  will  involve  AUSM  formulations  in  the  context  of  multi-dimensional  reconstruction. 

4.5  Ground  Based  Radiation  Data  (EAST  at  NASA  Ames) 

Johnston  [40]  presented  a  comparison  between  EAST  shock  tube  radiation  measurements  [61]  and  the  HARA 
radiation  model  [62,  63].  Equilibrium  and  nonequilibrium  radiation  measurements  were  studied  for  conditions 
relevant  to  lunar-return  shock-layers;  specifically  shock  velocities  ranging  from  9  to  1 1  km/s  at  initial  pressures 
of  0.1  and  0.3  Torr.  The  simulated  shock-tube  flow  was  assumed  one-dimensional.  The  shock-tube  flowfield 
was  simulated  using  the  LAURA  Navier-Stokes  solver  to  predict  the  flow  around  a  5  m  sphere  at  the  specified 
shock-tube  conditions.  A  two-temperature  thermochemical  nonequilibrium  11-species  (N,  N+,  NO,  NO+,  N2, 
N2+,  O,  0+,  O2,  02+,  and  e")  air  model  [64,  65]  was  applied  in  these  computations,  and  the  detailed  nonequi¬ 
librium  radiation  prediction  was  obtained  from  the  HARA  code  in  an  uncoupled  manner  (meaning  radiative 
coolong  was  not  accounted  for).  Similar  analyses  have  been  performed  using  other  flowfield  and  radiation 
solvers  by  Bose  et  al.  [66]  and  Panesi  et  al.  [67],  among  others. 

A  significant  fraction  of  the  radiative  heating  at  typical  peak-heating  lunar-return  conditions  is  emitted 
from  shock-layer  regions  in  thermochemical  equilibrium.  According  to  flowfield  predictions,  the  0.3  Torr 
EAST  cases  are  in  thermochemical  equilibrium  throughout  most  of  the  usable  test  time.  Thus,  they  provide  a 
valuable  means  for  validating  the  present  radiation  model  in  the  important  equilibrium  regime.  Comparisons 
between  experiment  and  predictions  in  this  regime  are  also  valuable  because  the  problems  of  nonequilibrium 
chemistry  and  non-Boltzmann  state-populations  do  not  complicate  the  theoretical  modeling,  therefore  allowing 
the  comparisons  to  focus  on  the  equilibrium  radiation  properties  and  not  the  nonequilibrium  rates. 

A  comparison  between  the  measured  and  predicted  intensity  spectrum  at  3.2  cm  behind  the  shock  is  shown 
in  Fig.  16(a)  for  the  10.34  km/s  case  at  0.3  Torr.  The  measured  wavelength  range  for  this  case  is  700  -  900 
nm,  with  the  dominant  spectral  features  being  the  atomic  line  multiplets.  The  significance  of  these  various 
spectral  features  is  indicated  in  this  figure  by  the  cumulative  wavelength-integrated  ( Jc)  curve,  which  has  been 
multiplied  by  100  for  scaling  purposes.  Note  that  a  constant  value  of  0.8  W/cm3/sr/^  was  subtracted  from  the 
measurements  for  this  case  to  force  the  continuum  radiation  to  agree,  as  discussed  in  [40].  This  subtracted 
value,  which  result  is  a  total  Jcxl00  of  14.4  W/cm3/sr,  is  interpreted  as  an  artifact  of  calibration  or  an  untreated 
radiation  mechanism  from  contamination  species.  The  Jc  predictions  at  shock- velocities  1.5%  below  the  spec¬ 
ified  value  are  shown  in  Fig.  16(a)  to  indicate  the  sensitivity  to  the  reported  uncertainty  in  the  measured  shock 
velocity.  The  agreement  between  the  measured  and  predicted  Jc  values  is  within  10%  throughout  this  wave¬ 
length  range.  Similar  agreement  is  seen  in  Fig.  16(b),  which  shows  Jc  values  from  a  sample  wavelength  range 
as  a  function  of  distance  behind  the  shock.  The  influence  of  the  1.5%  shock- velocity  uncertainty  is  shown  in 
these  figures  along  with  the  prediction  obtained  by  assuming  a  Boltzmann  distribution  of  electronic  states.  The 
Boltzmann  and  non-Boltzmann  predictions  are  seen  to  converge  in  the  regions  of  thermochemical  equilibrium, 
as  required  by  theory. 

The  prediction  of  nonequilibrium  radiation  is  dependent  upon  the  modeling  of  the  nonequilibrium  chemistry 
and  the  non-Boltzmann  population  of  the  radiating  states,  in  addition  to  the  oscillator  strengths  and  line-widths 
required  for  equilibrium  radiation  predictions.  The  nonequilibrium  chemistry  and  non-Boltzmann  models, 
which  are  assumed  to  be  uncoupled,  have  significantly  larger  uncertainties  than  the  equilibrium  radiation  prop- 
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(a)  Intensity  as  function  of  wavelength  and  cumulative  intensity  (b)  Integrated  intensity(800-830nm)  as  function  of  distance  behind 
3.2  cm  behind  shock  front.  shock. 

Figure  16:  Comparison  of  measured  and  computed  intensities  (700  -  900  nm)  behind  shock  front  in  EAST  facility  at  10.34  km/s 
and  0.3  torr.  Dashed  lines  show  predictions  if  specified  velocity  is  reduced  by  1.5%. 


(a)  Intensity  as  function  of  wavelength  and  cumulative  intensity  (b)  Integrated  intensity (850-8 80  nm)  as  function  of  distance  behind 
3.1  cm  behind  shock  front.  shock. 

Figure  17:  Comparison  of  measured  and  computed  intensities  (700  -  900  nm)  behind  shock  front  in  EAST  facility  at  9.165  km/s 
and  0.1  torr.  Dashed  lines  show  predictions  if  specified  velocity  is  increased  by  1.5%. 


erties.  These  large  uncertainties  are  due  mainly  to  a  lack  of  experimental  measurements  at  the  high-temperature 
conditions  of  interest  in  a  peak-radiative  heating  shock-layer.  The  EAST  measurements  therefore  provide  a 
valuable  set  of  data  for  assessing  the  non-Boltzmann  rate-model  assembled  by  Johnston  [63]  and  applied  in 
HARA.  These  measurements  may  also  be  used  to  validate  the  quasi-steady  state  assumption  [68]  applied  in 
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the  non-Boltzmann  models  solution  procedure,  which  allows  the  non-Boltzmann  model  to  remain  uncoupled 
from  the  flowfield  model,  as  opposed  to  more  complicated  coupled  approaches[69,  70].  Furthermore,  be¬ 
cause  the  electron  temperature  and  flowfield  number  densities  obtained  from  the  LAURA  flowfield  solver  drive 
the  non-Boltzmann  model,  these  nonequilibrium  data  provide  an  opportunity  to  validate  the  thermochemical 
nonequilibrium  models  applied  by  LAURA. 

The  study  of  the  0. 1  Torr  data  will  be  similar  to  the  study  of  the  0.3  Torr  data  presented  above.  However,  the 
presence  of  significant  regions  of  nonequilibrium  and  the  fact  the  present  cases  are  mostly  optically  thin  allow 
the  discussion  to  focus  on  the  nonequilibrium  radiation  (and  the  difference  between  the  Boltzmann  and  non- 
Boltzmann  models).  A  comparison  between  the  measured  and  predicted  intensity  spectrum  at  3. 1  cm  behind  the 
shock  is  shown  in  Fig.  17(a)  for  the  9.165  km/s  case  at  0.1  Torr.  The  measured  wavelength  range  for  this  case 
is  700  -  900  nm.  In  addition  to  the  EAST  data  and  present  non-Boltzmann  prediction,  this  figure  also  shows  the 
result  of  increasing  the  shock  velocity  by  1 .5%  or  by  assuming  a  Boltzmann  distribution  of  electronic  states.  It  is 
seen  that  the  baseline  non-Boltzmann  prediction  under-predicts  the  data  by  nearly  half.  The  significantly  better 
prediction  for  the  1.5%  increase  in  the  shock  velocity  suggests  that  the  1.5%  uncertainty  is  responsible  for  the 
some  underprediction  of  the  measurements.  The  differences  seen  between  the  non-Boltzmann  and  Boltzmann 
predictions  confirm  the  nonequilibrium  state  of  the  gas  at  this  location. 

Figure  17(b)  compares  the  spatial  variation  of  the  integrated  intensity  over  a  sample  spectral  range  for  the 
measurement  and  predictions  shown  in  Fig.  17(a).  The  non-Boltzmann  prediction  is  similar  in  shape  but  lower 
in  magnitude  than  the  measurements.  An  important  point  to  be  made  here  is  that  measurements  confirm  the  sup¬ 
pression  of  the  significant  Boltzmann-predicted  emission,  which  is  characteristic  of  the  non-Boltzmann  model 
in  regions  of  nonequilibrium.  Although  the  non-Boltzmann  predicted  nonequilibrium  radiation  is  larger  than 
the  equilibrium  radiation  (which  is  approached  near  z  =  4  cm)  for  these  cases,  the  contribution  of  nonequilib¬ 
rium  radiation  to  the  wall  radiative  flux  for  CEV  is  small  because  the  length  of  the  nonequilibrium  region  is  less 
than  2  cm  along  the  stagnation-line,  while  the  equilibrium  region  is  roughly  15  cm.  The  Boltzmann  predicted 
nonequilibrium  radiation,  on  the  other  hand,  is  large  enough  in  magnitude  to  contribute  significantly  to  the 
wall  radiative  flux,  even  with  the  very  small  length  of  the  nonequilibrium  region.  The  experimental  validation 
that  the  nonequilibrium  radiation  is  closer  to  the  non-Boltzmann  prediction  than  the  Boltzmann  prediction  is 
therefore  valuable,  because  it  confirms  that  the  nonequilibrium  radiation  contribution  is  much  less  than  that 
predicted  by  the  Boltzmann  model. 

4.6  Fire  II,  Coupled  Radiation  Only 

The  Fire  II  data  set  from  1967  remains  the  touchstone  for  validation  of  entry  heating  simulations  of  radiative 
heating.  [7 1-73]  The  vehicle  had  three  overlay ed  beryllium  heatshields.  A  fresh  surface  was  exposed  twice 
during  entry  after  the  outermost  heatshield  was  jettisoned  as  its  integrity  became  compromised  due  to  thermal 
load.  Consequently,  radiation  data  was  obtained  without  the  complications  of  ablation  products.  Metallic  sur¬ 
faces  are  generally  considered  to  be  strongly  catalytic  but  the  potential  formation  of  an  oxide  coating  introduces 
uncertainty  in  the  proper  formulation  of  this  boundary  condition.  Therefore  results  are  presented  using  bound¬ 
ing  approximations  for  a  super-catalytic  condition  (species  mass  fractions  on  the  surface  equal  to  free  stream 
values)  and  for  a  non-catalytic  condition  (zero  species  fraction  gradient  at  the  wall).  Surface  temperatures  were 
specified  according  to  flight  data. 

Instrument  locations  are  presented  in  Fig.  18.  The  radiometers  record  radiation  intensity  contribution  be¬ 
tween  0.4  and  6.2  eV  (window  limits).  The  calorimeters  record  total  convective  and  absorbed  radiative  heating. 
The  absorbed  radiative  energy  is  calculated  as  a  function  of  wavelength.  [74]  Simulations  at  two  trajectory  points 
are  intended  to  illustrate  capabilities  at  a  strong  thermochemical  nonequilibrium  condition  (1636  seconds)  and 
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Figure  18:  Sensors  on  Fire  II  vehicle,  circle  represents  radiometer  locations,  x  represents  calorimeter  locations. 


(c)  Radiation  intensity  distribution  at  1636  seconds. 


(d)  Radiation  intensity  distribution  at  1643  seconds. 


Figure  19:  Total  and  radiative  heating  distributions  around  Fire  II  vehicle  at  trajectory  points  1636  and  1643. 
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a  near-equilibrium  condition  (1643  seconds). 

The  computed  heating  rates  for  the  subject  points  are  compared  to  experimental  data  in  Fig.  19.  The 
total  heating  rates  measured  by  the  calorimeters  are  presented  in  Figs.  19  (a)  and  (b).  The  solid  black  line  in 
these  figures  indicate  the  uncoupled  total  heating  (radiative  heating  computed  from  converged,  non-radiating 
case)  assuming  a  super-catalytic  boundary  condition.  The  solid  red  line  indicates  the  total  heating  with  full 
coupling  of  radiative  energy  transfer.  The  coupling  reduces  the  total  heating  by  5  to  10  %  for  these  trajectory 
points.  The  dashed  lines  indicate  the  coupled,  non-catalytic  heating  rate.  The  catalytic  effect  is  most  prominent 
(reduction  by  1/3)  at  the  strong  nonequilibrium  point  at  1636  seconds.  The  catalytic  effect  is  still  evident  at 
1643  seconds  with  a  net  reduction  of  5  %.  The  experimental  data,  indicated  by  blue  circles,  is  fully  bounded 
by  the  extremes  in  modeling  of  surface  catalysis.  The  uncoupled  radiation  intensity  is  presented  with  a  dashed 
line  and  the  coupled  radiation  intensity  is  presented  with  a  solid  red  line  in  Figs.  19  (c)  and  (d).  The  uncertainty 
in  the  experimental  data  in  these  figures  (blue  circles  and  error  bar)  bounds  the  computed  coupled  results.  The 
simulations  indicate  that  radiation  coupling  is  required  to  match  experimental  data  and  that  the  suite  of  physical 
models  used  to  compute  these  cases  are  consistent  with  experimental  data. 

4.7  Mars  Return,  Coupled  Radiation  and  Ablation 


Figure  20:  Temperature  profiles  for  the  Mars  return  case  with  and  without  coupled  radiation. 


A  benchmark  case  for  modeling  strongly  radiating  and  ablating  flowfields  consists  of  a  3.05  m  radius  sphere 
at  Voo  -  15.24  km/s  and  =  0.000255  kg/m3.  This  case,  which  represents  Mars  return  to  Earth,  has  been 
extensively  studied  for  over  30  years,  beginning  with  the  then  state-of-the-art  viscous  shock  layer  (VSL)  code 
methodology  [75-77].  In  Ref.  [78],  these  previous  studies  were  used  to  assess  the  recently  implemented 
coupled  radiation,  coupled  ablation,  and  free-energy  minimization  capability  of  the  LAURA  code  [35].  The 
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free-energy  minimization  capability  was  applied  to  provide  chemical  equilibrium  flowfield  solutions  consistent 
with  previous  VSL  solutions.  The  ablation  rates  and  chemical  composition  assumed  in  the  previous  studies 
were  applied,  and  the  resulting  stagnation  line  radiation,  temperature,  and  elemental  composition  were  shown 
to  agree  with  the  past  studies. 


Figure  21 :  Radiative  flux  to  the  wall  predicted  with  and  without  coupled  radiation  for  the  Mars  return  case. 


Figure  22:  Mass  fraction  of  two  dominant  ablation  products,  CO  and  C,  for  the  Mars  return  case  with  coupled  ablation  and 
radiation. 


Since  the  publication  of  Ref.  [78],  three  primary  advancements  have  been  made  to  the  study  of  this  Mars 
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return  example.  The  first  of  these  advancements  is  that  the  assumption  of  chemical  equilibrium  has  been  re¬ 
moved,  and  instead  the  modern  two-temperature  thermochemical  nonequilibrium  model  is  applied.  The  second 
of  these  advancements  is  that  the  ablation  rate  is  computed  assuming  steady-state  ablation  as  part  of  the  flow- 
field  solution,  and  therefore  the  specification  of  the  ablation  rate  is  not  required.  The  procedure  required  for 
this  ablation  rate  computation  is  described  by  Johnston  et  al.  [79].  Assuming  the  elemental  mass  fractions  of 
injected  gas  are  [C,  H,  O,  N]  =  (0.92,  0.022,  0.049,  0.009),  the  predicted  ablation  rate  at  the  stagnation  point 
for  the  present  Mars  return  case  is  0.445  kg/s/m2.  The  third  of  these  advancements  is  the  treatment  of  the 
sphere’s  afterbody,  which  allows  for  the  treatment  of  the  wake.  For  the  present  results,  the  surface  of  the  sphere 
afterbody  was  assumed  to  be  non-ablating,  radiative  equilibrium,  and  fully-catalytic  to  homogenous  recombi¬ 
nation.  This  assumption  of  a  non-ablating  afterbody  allows  the  influence  of  forebody  ablation  products  on  the 
afterbody  to  be  studied. 

The  influence  of  coupled  radiation,  without  coupled  ablation,  on  the  vibrational-electronic  temperature  at 
various  points  in  the  flowfield  is  shown  in  Fig.  20.  The  radiative  heating  at  the  surface  for  these  cases  is  shown 
in  Fig.  21.  These  figures  show  that  the  radiatively  cooled  gas  from  the  strongly  radiating  stagnation  region 
flows  downstream  and  reduces  the  temperatures  in  the  weakly  radiating  downstream  regions  of  the  flow.  This 
non-local  radiative  cooling  effect  causes  correlations  such  as  the  Tauber- Wakefield  approximation[80],  which 
is  dependent  upon  the  local  radiative  flux,  to  under-predict  the  radiative  cooling  effect  in  dowstream  regions. 
This  effect  is  apparent  in  Fig.  21,  especially  in  the  wake  where  the  Tauber- Wakefield  approximation  predicts 
very  little  cooling. 

The  mass  fractions  of  CO  and  C  predicted  throughout  the  flowfield  by  assuming  an  ablating  forebody  and 
non-abalting  afterbody  (as  mentioned  previously)  are  shown  in  Fig.  22.  This  figure  shows  that  although  both 
CO  and  C  are  formed  as  a  result  of  the  forebody  ablation  products,  they  flow  into  the  low  pressure  wake  and 
dominate  the  species  composition  near  the  body.  The  influence  of  ablation  on  the  radiative  heating  is  shown 
in  Fig.  23.  On  the  forebody,  a  reduction  in  the  radiative  heating  similar  to  that  shown  previously  in  Ref.  [78] 
is  seen,  while  the  forebody  shows  a  surprising  increase  in  the  radiation  with  the  introduction  of  ablation.  This 
increase  is  a  result  of  emission  from  the  CO  4+  band  system,  whose  contribution  to  the  radiative  flux  in  the 
wake  is  non-negligible  because  of  the  large  CO  concentration  in  that  region. 


s(m) 


s(m) 


Figure  23:  Radiative  flux  to  the  wall  predicted  with  and  without  coupled  ablation  for  the  Mars  return  case. 
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4.8  Apollo  4,  Coupled  Radiation  and  Ablation 

The  Apollo  4  test  case,  recently  reviewed  by  Park[8 1],  provides  flight  data  on  radiative  heating  including  effects 
of  ablation  products.  Free  stream  conditions  are  =  10.25  km/s  and  =  0.000341  kg/m3  on  a  3  meter  radius 
sphere.  Elemental  mass  fractions  of  injected  gas  are  [C,  H,  O,  N]  =  (0.63,  0.06,  0.30,  0.01).  These  fractions 
were  derived  from  data  discussed  by  Park[81]  in  which  quasi-steady  ablation  is  assumed  and  the  elemental 
fractions  of  char  and  pyrolysis  gas  are  combined  according  to  the  fraction  of  their  individual  blowing  rates. 
The  wall  temperature  is  set  to  2500  K  and  the  dimensionless  blowing  rate  m  is  0.0086  [81].  Thermochemical 
nonequilibrium  model  is  applied.  The  case  was  initialized  from  a  non-ablating  solution  in  which  the  equilibrium 
catalytic  wall  boundary  condition  was  applied. 


(a)  Intensity  integrated  over  frequency  over  two  paths 


Figure  24:  Comparisons  with  Apollo  4  radiometer  data. 


Fig.  24  (a)  plots  the  cumulative  intensity  arriving  at  the  radiometer  starting  from  0.4  eV  and  concluding 
at  6.2  eV  (the  radiometer  window  is  opaque  outside  this  range)  .  The  path  length  indicated  by  the  solid  line 
includes  the  shock-layer  thickness  immediately  above  the  radiometer.  The  path  length  indicated  by  the  dashed 
line  includes  the  shock-layer  thickness  plus  8  cm  of  cavity  depth  (Fig.  24  (b))  which  is  assumed  to  be  filled 
with  gas  at  conditions  on  the  surface.  (There  was  no  flow  field  simulation  that  included  a  cavity.)  From  0.4  to 
approximately  2.0  eV  there  is  essentially  no  difference  in  the  intensity  reaching  the  radiometer  over  either  path 
length.  The  intensity  reaching  the  radiometer  above  2eV  becomes  path  length  dependent,  indicating  that  some 
of  the  radiation  above  2eV  is  being  absorbed  by  the  gas  in  the  cavity.  The  circular  symbol  with  uncertainty  bar 
indicates  the  measured  cumulative  intensity.  It  is  convenient  to  put  this  symbol  at  the  end  of  the  window  range 
at  6.2  eV  but  the  data  in  fact  represents  the  integrated  intensity  from  0.4  to  6.2eV.  This  figure  shows  that  for 
the  assumed  blowing  conditions  the  predicted  cumulative  intensity  at  the  surface  exceeds  the  measured  value 
by  approximately  20  %.  However,  if  the  effect  of  absorption  over  a  longer  path  length  through  the  cavity  is 
included  then  excellent  agreement  with  measurement  is  obtained. 
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Figure  25:  Apollo  4  flowfield  with  coupled  radiation  and  ablation. 


distance  along  centerline  (m) 


Afterbody 


Figure  26:  Influence  of  coupled  ablation  and  radiation  on  the  forebody  and  afterbody  convective  heating  for  the  Apollo  4  case. 


We  caution  that  this  agreement  does  not  constitute  validation.  The  specified  blowing  rates  and  elemental 
mass  fraction  of  the  ablation  products  are  estimates  [81]  -  not  computed  on  any  first  principles  within  this 
simulation.  Still,  the  results  confirm  that  for  a  reasonable  estimate  of  boundary  conditions  there  is  consistency 
with  measured  data  -  providing  more  supporting  evidence  that  the  coupled  radiation  and  ablation  algorithms 
are  correctly  implemented.  Further  details  of  this  simulation  with  more  consideration  of  parametric  variation  is 
available  in  the  literature.  [78] 

Unlike  the  study  of  Ref.  [78]  discussed  above,  which  approximated  the  Apollo  stagnation  region  with  a 
3  m  hemisphere,  the  present  results  are  obtained  using  the  actual  Apollo  geometry  at  angle  of  attack.  Ele- 
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Figure  27:  Influence  of  coupled  ablation  and  radiation  on  the  forebody  radiative  heating  for  the  Apollo  4  case. 


mental  mass  fractions  of  the  char  and  pyrolysis  gas  are  taken  from  Bartlett  et  al.  [82]  and  assume  mechanical 
silica  removal.  Steady-state  ablation  is  assumed  for  the  computation  of  the  pyrolysis  injection  rates  and  wall 
temperature.  (Recall  in  the  previous  case  the  surface  temperature  and  ablation  rate  was  specified.)  This  as¬ 
sumption  is  appropriate  for  the  peak  heating  conditions  considered  here.  The  flowfield  was  computed  assuming 
two-temperature  thermochemical  nonequilibrium.  Similarly  to  the  Mars-return  case,  only  the  forebody  is  as¬ 
sumed  to  be  ablating.  The  afterbody  is  assumed  to  be  fully-catalytic  for  homogenous  recombination  and  the 
temperature  is  obtained  from  the  assumption  of  radiative  equilibrium. 

The  vibrational-electronic  temperature  throughout  the  flowfield  is  shown  in  Fig.  25  (a).  This  flowfield 
includes  coupled  radiation  and  ablation.  The  mass  fractions  of  the  dominant  ablation  product,  CO,  are  shown  in 
Fig.  25  (b).  Note  that  the  CO,  which  originates  on  the  forebody,  flows  into  the  afterbody  and  result  in  up  to  30% 
CO  on  the  afterbody  surface.  The  influence  of  this  ablation  product  contamination  is  shown  in  Fig.  26  to  result 
in  a  slight  decrease  in  the  convective  heating  on  the  wind- side  afterbody.  For  lower  velocity  entry  conditions, 
this  decrease  has  been  found  to  be  as  large  as  20%.  The  decrease  in  the  convective  heating  seen  on  the  forebody 
in  Fig.  26  is  significant  because  this  region  is  locally  ablating.  The  influence  of  coupled  radiation  is  shown  in 
this  figure  to  have  a  small  impact  on  the  convective  heating.  The  influence  of  coupled  ablation  and  radiation 
on  the  forebody  radiative  heating  (the  afterbody  radiation  is  negligible)  is  shown  in  Fig  27.  The  influence  of 
coupled  radiation  is  shown  to  reduce  the  radiative  heating  by  nearly  30%,  while  the  influence  of  ablation  reduces 
it  by  only  about  5%.  Note  that  the  ablation  rates  predicted  for  this  case  are  an  order-of-magnitude  less  than  for 
the  Mars  return  case  discussed  previously,  in  which  case  ablation  had  a  significant  impact  on  the  radiation. 

4.9  Towed  Ballute 

A  series  of  wind  tunnel  tests  were  conducted  at  Langley  in  Mach  6  air  and  CF4  to  gather  surface  heating  data 
(thermophosphor)  and  Schlieren  video  to  observe  unsteady  motion.  [16]  Some  of  these  tests  were  exploratory  in 
nature,  examining  ways  in  which  flexible,  towed  decelerators  behaved  dynamically  in  a  ground  based  facility  - 
accepting  the  fact  that  Mach  number  and  dynamic  pressure  are  far  from  the  flight  condition.  We  discuss  here 
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(a)  Run  18  -  Schlieren 

Figure  28:  Single  frame  of  experiment  and  CFD  simulation  for  Run  18  with 


-4  -2  0  2 


z 

(b)  Simulation  of  Run  18  -  density  contours 

weakly  unsteady  flow  (small  amplitude  excursions). 


a  test  of  a  hemispherically  capped  cylinder  extending  through  the  center  of  a  toroid.  The  spacecraft  (cylinder) 
diameter  is  0.50  inches.  The  toroid  diameter  is  6  inches.  It  produces  a  bow  shock  which  is  directed  just  inside 
the  inner  boundary  of  the  toroid.  As  the  cylinder  is  pushed  further  upstream,  unsteady  motion  ensues.  As  the 
cylinder  is  retracted  the  motion  appears  steady  (at  least  with  regard  to  images  captured  in  the  Schlieren).  The 
topology  of  the  wind  tunnel  case  is  different  than  the  flight  system  in  that  the  cylindrical  sting  is  not  present  in 
flight.  Nevertheless,  the  same  behavior  of  the  recirculating  flow  moving  forward  due  to  the  compression  in  the 
core  of  the  toroid  is  inferred  from  the  pulsing  bow  shock  behind  the  nose  of  the  cylinder. 

Figure  28(a)  is  a  single  frame  from  a  Schlieren  video  of  Run  18  in  the  20-Inch  Mach  6  Air  Tunnel  at  a 
Reynolds  number  of  1  million  per  foot.  The  distance  from  the  center  of  the  hemispherical  nose  to  the  midplane 
of  the  toroid  is  3  inches.  The  video  frame  rate  is  not  sufficiently  high  to  resolve  a  single  cycle  of  recirculating 
flow.  Instead,  it  captures  instantaneous  (to  exposure  time)  locations  of  the  bow  shock  as  it  responds  to  the 
unsteady,  separated  flow  beneath  it.  The  observed  range  of  motion  in  this  case  is  characterised  as  incipient 
unsteady  response.  There  is  a  slight  fluttering  with  small  wavelength  visible  in  the  bow  shock  image.  This 
incipient  case  is  considered  an  important  test  because  we  want  to  confirm  the  ability  of  LAURA  to  identify  the 
presence  and  extent  of  unsteady  motion. 

A  single  frame  from  the  simulation  of  Run  18  with  LAURA  is  presented  in  Fig.  28(b)  showing  density 
contours.  In  the  movie,  there  is  significant  movement  of  the  reflected  shock  off  the  cylinder  which  drives  the 
recirculation  flow  back  to  the  nose.  These  structures  cannot  be  picked  up  in  the  Schlieren.  The  response  of  the 
bow  shock  over  the  cylinder  is  more  subdued.  The  fluttering  of  the  bow  shock  in  the  CFD  simulation  appears 
to  be  somewhat  stronger  than  observed  in  the  experiment.  The  key  point  here  is  that  CFD  in  fact  picked  up  an 
unsteady  response  when  one  is  observed  in  the  experiment.  In  contrast,  Run  16  (not  shown  here)  in  which  the 
cylinder  was  retracted  another  inch  backward  produced  a  flow  which  was  steady  in  the  Schlieren  and  steady  in 
the  CFD  simulation.  This  agreement  builds  credibility  for  the  first  time  that  the  available  tools  are  at  least  able 
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(a)  Run  20  -  Schlieren 

Figure  29:  Single  frame  of  experiment  and  CFD  simulation  for  Run  20  with  strong  unsteady  flow  (large  amplitude  excursions). 
The  magnitude  of  oscillations  is  qualitatively  captured. 


to  discern  unsteady  motion  when  it  occurs.  The  more  complex  task  of  validating  the  time  accurate  response  is 
of  less  concern  in  that  we  want  to  avoid  the  unsteady  motion  in  the  flight  system  in  the  first  place  (or  at  least 
insure  that  the  extent  of  unsteady  motion  does  not  compromise  the  flight  system). 

An  example  of  a  violently  unsteady  response  which  would  compromise  system  integrity  is  evident  in  Run 
20  in  which  the  cylinder  is  pushed  an  additional  inch  forward  (as  compared  to  Run  18).  In  this  configuration 
(Fig.  29(a))  the  bow  shock  impinges  on  the  toroid.  A  more  complex  wave  pattern  reflects  off  the  cylinder  in  the 
core  of  the  toroid,  and  a  massively  unsteady  response  ensues.  The  bow  shock  off  the  cylindrical  nose  moves 
from  the  inner  to  outer  edge  of  the  toroid  -  a  situation  totally  unacceptable  in  flight  from  both  a  heating  and 
system  stability  perspective.  The  unsteady  simulation  of  Run  20  (Fig.  4.9)  shows  an  equivalent  amplitude  of 
motion  as  the  experiment.  The  bow  shock  off  the  spherically  capped  cylinder  swings  from  the  inside  to  outside 
edge  of  the  trailing  toroid. 

5.0  UNCERTAINTY  ANALYSES 

Numerical  simulations  should  be  expected  to  publish  uncertainties  in  the  same  way  uncertainties  in  experi¬ 
mental  data  are  required.  Simulations  generally  have  two  sources  of  error:  (1)  truncation  error  associated  with 
spatial,  temporal,  and  spectral  discretization  and  (2)  modeling  error  associated  with  turbulence,  kinetics,  trans¬ 
port  properties,  thermal  relaxation,  radiation,  and  boundary  conditions.  Truncation  error  may  be  quantified 
through  grid  convergence  studies.  Application  of  output  based  grid  convergence  (adjoint  based  adaptation)  has 
been  used  to  drive  grid  adaptation  to  meet  user  specified  uncertainties  in  output  quantities.  Modeling  error 
may  be  characterized  as  model  form  uncertainty  and  model  parameter  uncertainty.  Quantification  of  model 
form  uncertainty  (e.g.  algebraic  turbulence  model  vs  one-  or  two  equation  turbulence  model)  requires  extensive 
experimental  validation  and  the  magnitudes  are  usually  dependent  on  local  flow  conditions  (mission  specific, 
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application  specific).  Quantification  of  uncertainty  derived  from  model  parameter  uncertainties  in  computa¬ 
tional  aerothermodynamic  simulations  has  been  discussed  by  Kleb  et.  al.  [83,  84].  In  [83]  a  markup  format  to 
input  data  files  is  described  to  automate  uncertainty  analyses  in  a  simulation.  (See  Fig.  30)  In  [84]  estimates  of 
epistemic  uncertainties  in  radiation  model  parameters  to  HARA  are  defined  and  the  net  uncertainty  to  several 
test  problems  are  calculated  —  including  the  FIRE  II  simulation  presented  previously.  The  intensity  between  0 
and  6  eV,  above  200  nm,  was  measured  by  the  Fire  II  radiometer  and  was  presented  in  Fig.  19.  The  uncertainty 
analysis  results  in  intensity  of  4.8  ±  1.5  W/cm2/sr  for  the  1636  case,  which  compares  well  with  the  measured 
result  of  4.9  ±1.1  W/cm2/sr.  For  the  1643  case,  the  analysis  predicts  51.1  ±  13.8  W/cm2/sr,  which  overlaps 
the  error  bounds  of  the  measured  result  of  63.0  ±  14.5  W/cm2/sr. 


- kinetics_data - 

- kinetics_data  fzy - 

02  +  M  <=>  20  +  M 

02  +  M  <=>  20  +  M 

1.000e+22  -1.50  5.950e+04 

1.000e+22  -1.50  5.950e+04 

teffl  =  2 

teffl  =  2 

expl  =  0.7 

expl  =  0.7+/-0.25U  'Eta_02' 

t  eff  min  =  1000. 

t  eff  min  =  1000. 

t  eff  max  =50000. 

—  markup  — ► 

t_eff_max  =  50000. 

02  =  0.2 

02  =  0.2+/-1O  ' A_02+02 ' 

N2  =  0.2 

N2  =  0.2+/-lo  ' A_02+N2 ' 

NO  =  0.2 

NO  =  0.2+/-lo  ' A_02+NO ' 

N  =1.0 

N  =  1.0+/-1O  ' A_02+N ' 

0  =  1.0 

0  =  1.0+/-1O  ' A_02+0 ' 

Figure  30:  Formats  for  marking  up  parameter  uncertainties  on  physical  property  input  files  to  LAURA. 


6.0  CONCLUDING  REMARKS 

Mission  requirements  for  entry,  descent,  and  landing  (EDL)  of  high  mass  at  Mars  and  exploration  of  the  outer 
planets  highlight  the  need  for  new,  innovative  vehicle  designs.  Depending  on  the  mission,  simulations  may 
require  coupling  of  some  combination  of  radiation,  massive  ablation,  surface  deflections,  and  jet  interactions 
-  sometimes  in  a  time-accurate  manner.  Addressing  validation  of  the  requisite  suite  of  physical  models  in 
computational  aerothermodynamics  (CA)  remains  a  high  priority.  The  need  to  explore  innovative  solutions  to 
design  problems  has  become  an  equally  high  priority.  The  design  space  for  EDL  systems  (e.g.  layout  and 
thrust  schedules  of  retro-rockets,  attachment  details  for  trailing  ballutes,  effect  of  deployable  infrastructure  on 
localized  heating  and  unsteady  interactions,  effects  of  large  shape  change  on  lifting  bodies)  require  that  CA  be 
brought  into  the  design  process  early.  These  design  challenges  tend  to  introduce  highly  non-linear  interactions 
that  are  not  amenable  to  simple  engineering  analysis  tools.  Challenges  within  CA  discussed  herein  have  been 
framed  in  the  context  of  simulation  for  innovation.  Simulation  for  innovation  has  three  defining  features:  (1) 
conceptual  systems  can  be  quickly  set  up  for  simulation  and  optimization,  (2)  robust  grid  adaptation  to  a  user 
specified  discretization  tolerance  is  available,  and  (3)  credible  estimates  of  simulated  aerothermal  environment 
uncertainties  are  published.  Obstacles  to  any  of  these  features  impedes  innovative  design  -  be  it  an  EDL  system 
to  land  tens  of  metric  tons  on  the  Martian  surface  or  set-up  of  an  experiment  for  validation  in  a  ground-based 
facility. 

Simulation  for  innovation  in  CA  is  the  goal.  Two  algorithmic  topics  that  must  be  addressed  to  achieve  this 
goal  are  discussed.  The  first  topic  addresses  the  need  to  obtain  accurate  heating  simulation  on  grid  systems  that 
are  easily  adapted  to  complex  flow  topologies.  Tetrahedral  grids  provide  the  greatest  flexibility  for  adaptation 
but  are  poorly  suited  for  prediction  of  heating  and  shear  in  hypersonic  flows  using  conventional  finite- volume 
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algorithms.  A  recent  development  in  multi-dimensional,  inviscid  flux  reconstruction  substantially  overcomes 
these  limitations  as  demonstrated  in  two  simple  challenge  problems  as  well  as  two  more  complex  configurations 
—  double  cone  and  shuttle  orbiter.  The  second  topic  addresses  the  infrastructure  needs  to  support  fully  coupled 
ablation,  radiation,  and  shape  change  in  a  CA  simulation.  A  free  energy  minimization  option  enables  equilib¬ 
rium  catalytic  wall  and  equilibrium  ablation  boundary  conditions  as  well  as  equilibrium  shock  layer  simulations 
in  which  elemental  mass  fractions  vary  due  to  ablation.  The  HARA  radiation  code,  using  the  smeared-rotational 
band  (SRB)  method  for  molecular  band  radiation,  requires  only  30  to  50  seconds  compute  time  per  line  of  sight. 
This  relative  fast  computation  time  enables  coupled  radiation  with  the  flow  solver  for  3D  configurations. 

Examples  of  code  assessment  and  validation  using  the  structured  flow  solver  LAURA  and  the  unstructured 
flow  solver  FUN3D  in  problems  of  particular  interest  to  the  CA  community  are  presented.  These  cases  illustrate 
important  simulation  capabilities  but  also  serve  to  illustrate  where  progress  is  needed  to  achieve  simulation 
for  innovation.  Both  ground-based  and  flight  data  are  computed.  Simulations  of  shuttle  orbiter  heating  on 
both  the  structured,  cell-centered  LAURA  code  formulation  and  on  the  unstructured,  node-based  FUN3D  code 
formulation  include  important  effects  of  surface  catalysis.  Simulations  of  radiative  energy  transfer  at  near¬ 
equilibrium  and  non-equilibrium  conditions  compare  well  with  measured  flight  data  (FIRE  II)  and  ground- 
based  shock  tube  data  in  EAST.  Coupled  massive  ablation  ( m  —  0.2)  and  radiation  are  tested  on  a  classic  Mars 
return  problem  with  good  comparisons  to  previous  solutions  using  viscous  shock  layer  methods.  Comparison  to 
radiometer  data  on  Apollo  4  test  flight  including  coupled  ablation  and  radiation  shows  good  agreement  though 
the  blowing  rate  in  this  case  is  prescribed.  Simulation  of  a  ground-based  test  of  supersonic  retro-propulsion 
demonstrates  the  difficulty  and  criticality  of  grid  adaptation  as  judged  by  the  plume  shape  and  terminal  shock 
shape.  Simulation  of  unsteady  flow  as  a  function  of  relative  position  of  the  spacecraft  and  towed  ballute  shows 
qualitative  agreement  with  ground-based  tests.  In  every  case,  significant  expert  user  intervention  was  required 
to  craft  a  grid  and/or  to  tune  the  infrastructure  for  multi-physics  coupling. 

Finally,  in  an  initial  effort  to  automate  the  publishing  of  uncertainties  in  a  given  simulation,  an  algorithm  is 
briefly  described  that  allows  one  to  compute  the  net  uncertainty  on  output  quantities  as  a  function  of  epistemic 
uncertainties  in  physical  model  parameters.  An  example  involving  radiative  heating  is  provided.  Though  this 
algorithm  is  unable  to  formally  address  model  form  uncertainties  it  is  considered  an  important  step  in  realizing 
the  third  element  of  simulation  for  innovation  —  credible  estimates  of  simulated  aerothermal  environment 
uncertainties  are  published. 

7.0  APPENDIX  A  -  MULTI-DIMENSIONAL  RECONSTRUCTION 

The  multi-dimensional  reconstruction  algorithm  used  to  enable  a  tetrahedral  grid  system  for  surface  heating  is 
provided  here  for  the  convenience  of  the  reader.  More  details  are  available  in  the  original  documentation.  [31] 
The  new  algorithm  is  implemented  within  the  flow  solver  FUN3D  -  a  node  based,  fully  unstructured,  finite- 
volume  solver  of  the  Euler  and  Navier-Stokes  equations. [85]  The  method  uses  a  Green-Gauss  formulation  of 
gradients  at  the  nodes  for  multi-dimensional,  inviscid  flux  reconstruction.  Viscous  gradients  are  also  computed 
from  a  Green-Gauss  formulation.  A  suite  of  modules  includes  all  of  the  gas  physics  models  in  LAURA  and 
VULCAN[86]  for  thermodynamics,  transport  properties,  chemical  kinetics,  and  thermal  relaxation.  The  base¬ 
line  inviscid  flux  reconstruction  algorithms  utilize  quasi-one-dimensional  (edge-based)  reconstruction  using 
Roe’s  averaging  with  Yee’s  symmetric  total  variation  diminishing  algorithm  adapted  for  unstructured  grids. 

A  triangular  grid  system  is  presented  in  Fig.  31.  The  two-dimensional  reconstruction  algorithm  is  developed 
relative  to  this  figure.  The  extension  of  the  algorithm  to  three  dimensions  is  presented  in  Appendix  A.  This 
derivation  is  repeated  from  the  original  source  [29]  for  the  convenience  of  the  reader  with  special  emphasis 
placed  on  updates  to  the  algorithm. 
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Assume  a  node  based  scheme  as  illustrated  in  Fig.  31  with  (x,  y)  coordinates  of  all  nodes  (1,2,3, . . . )  given 
and  dependent  variables  q  at  these  same  nodes  to  be  determined  by  solution  of  conservation  laws.  The  conser¬ 
vation  laws  are  discretized  relative  to  the  dual  control  volume,  represented  by  dashed  lines  in  Fig.  31.  The  dual 
volume  around  node  3  passes  through  the  centroids  of  all  elements  surrounding  the  node  (A,  F>,  C,  E ,  F ) 
and  the  midpoints  of  all  edges  separating  these  elements. 


Figure  31 :  Schematic  of  node-based,  unstructured  grid  system  showing  element  centroids  (uppercase  letters),  edge  midpoints 
(lowercase  letters),  nodes  (numerals),  the  dual  volume  (dashed  lines),  and  the  principal  directions  (dotted  lines  with  arrowhead). 


7.0.1  Reconstruction  Overview  -  Principal  (Rotated)  Direction 

Consider  a  conventional,  quasi-one-dimensional,  first-order  reconstruction  of  inviscid  flux  across  face  AcB 
separating  the  dual  volumes  around  nodes  1  and  3.  This  flux  is  defined  as  a  function  of  information  at  nodes  1 
and  3  of  edge  c: 

f AcB  =  1  [f3  +  fi  -  Rc|Ac|R“1(q3  -  qi)]  (20) 

Inviscid  fluxes  are  then  defined  using  a  loop  over  edges  and  the  reconstruction  direction  is  fully  determined  by 
the  grid. 

The  basic  idea  in  multi-dimensional  reconstruction  is  that:  (1)  reconstruction  directions  should  be  based  on 
local  flow  characteristics  and  not  be  constrained  by  the  grid;  and  (2)  the  flux  components  for  any  face  will  utilize 
a  stencil  employing  all  nodes  of  a  surrounding  element.  Conventional,  node-based  schemes  loop  over  edges 
when  using  quasi-one-dimensional  reconstruction.  In  contrast,  the  new  multi-dimensional  reconstruction  loops 
over  elements.  Consider  the  two-dimensional  element  defined  by  nodes  (1,  2, 3)  with  centroid  A  in  Fig.  31. 
The  computation  of  inviscid  flux  in  this  element  employs  the  same  infrastructure  as  the  viscous  flux. 

A  principal  direction  in  the  element,  xf ,  is  defined  by  the  direction  of  Vp,  computed  using  the  same  Green 
Gauss  evaluation  of  the  derivatives  that  are  already  used  to  compute  the  viscous  terms.  If  density  is  locally 
constant  then  the  principal  direction  aligns  with  the  local  velocity  computed  as  the  average  of  nodal  values. 
Inviscid  flux  within  the  element  will  be  computed  along  and  orthogonal  to  the  principal  flow  direction. 

The  flux  components  are  defined  using  nodal  weights  based  on  the  Green-Gauss  formulation  for  derivatives 
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in  the  component  directions.  Thus, 


dU 

dx' 


1 

£Ia 


k= faces 


(21) 


where  is  the  average  value  of  f  across  all  nodes  defining  surface  k  of  element  A.  With  reference  to  Fig.  31 
this  derivative  can  be  written 


dfx',A 

dx' 


a3,x'(h,x'  +  fl,x')  +  al,x'(^3,x'  +  h,x')  +  a2,x'(h,x'  +  ^3,x') 


(22) 


where 


^k,x' 


^k^k^x' 
2  04 


and  k  is  a  face  index  identified  by  the  number  of  the  opposite  node.  The  sum  J2k= faces  ak,x' 
closed  element.  It  is  convenient  to  reorder  and  scale  the  coefficients  as  follows. 

(a2,x'  +  a3,x')fl  +  (a3,x'  +  al,x')h  +  +  OL2,x  0^3 

\Pl,x'fl  +  P 2,x'h  +  /Vf3] 

Note  that  Ax'  can  be  interpreted  as  a  distance  across  the  element  in  x'  direction.  It  is  defined 

^^7  =  (\an-l,x'  +  an+l,x'  |) 

n=nodes 


dx' 


(23) 
0  for  any 


(24) 

(25) 


(26) 


where  a  cyclic  indexing  is  assumed.  Note  that  the  original  formulation  [29]  is  corrected  here  with  a  factor 
2  in  the  numerator  which  provides  the  distance  between  virtual  nodes.  The  reordering  and  scaling  yield  the 
following  relations  for  Prux> . 


Pn,x' 

Ax  ( OLn—\  x’  T"  OLn- fl^') 

(27) 

^  ^  Pn,x' 

Tl 

=  0 

(28) 

1  Pn,x'  | 

=  2 

(29) 

n 


Two  options  for  computing  ix>  are  developed.  Both  mirror  the  baseline  quasi-one-dimension  reconstruction 
using  Roe’s  averaging  across  nodes.  The  first  option  averages  dependent  variables  at  two  virtual  nodes  used  in 
the  reconstruction.  The  second  option  may  be  considered  a  weighted  average  of  a  rotated  difference  scheme  on 
existing  edges.  No  virtual  nodes  are  created. 


7.0.2  Option  1  -  Virtual  Node  Averaging 

A  right  and  left  virtual  state  for  computing  are  computed  as  follows. 


*T R,x'  ^  ^  (l/^n^'l  “1“  Pn,x')Q.n 

(30) 

n=nodes 

*\L,x'  ~  ^  ^  (IAt,,#'!  —  Pn,x')^Ln 

(31) 

n=nodes 
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l  r 


f^'  =  - 


2  L 


\-^-Um,x'\(dc{xi  (iq; 


ix'  ,lim) 


(32) 


where  the  left  and  right  states  of  f  are  functions  of  the  left  and  right  states  of  q,  respectively.  The  matrices 
Hxf  and  Aum,x'  are  computed  as  functions  of  the  Roe’s  average  of  the  right  and  left  states.  (Forming  q#^/ 
and  q.L,x'  from  averages  of  1  / u,v,w ,  1/T  have  provided  the  most  robust  simulations  as  judged  by  level  of 
temperature  undershoots  across  strong  bow  shocks.) 


dqx'  —  -F ^x'{qR,x'  O.L,x') 

^  ^  fin,x'Qjri 
n=nodes 


(33) 


dq.R,x'  —  ^  ^  dAi,#'!  H"  Pn,x'  Qn,GG 

n=nodes 

dc[L,x'  =  Ax  R^/  ^  ^  (l/^n^'l  fin,x'  Q.n,GG 


n=nodes 

dqx',lim  =  $  rninmod 
where  <f>  is  an  auxiliary  limiter  (0  <  $  <  1)  defined  in  a  later  section. 


^dc[R  xf ,  2dqx/,  2dq^  ^  {dqRx/  +  dqRx 


(34) 

(35) 

(36) 


7.0.3  Option  2  -  Weighted  Average  of  Edges  to  Principal  Node 


Select  the  largest  of  |  /3njXf  |  to  identify  the  principal  node  for  construction  of  fx>  and  express  it  as  a  function  of 
the  remaining  coefficients.  Assume  node  3  is  the  principal  node.  Then 

dfx'  =  (3i,x'fi  +  fay  f 2  +  fi^x'h 

—  Pl,x'fl  +  P2,x'f2  —  (Pl,x'  +  ^2,x')f3  (37) 

=  Pl,x'{fl  ~  f?>)  +  P2,x'{f2  ~  h) 

The  reconstructed  flux  in  the  direction  x ix>,  is  computed  as  a  weighted  average  of  surrounding  edges.  Fur¬ 
thermore,  it  is  noted  that  if  any  of  the  surrounding  edges  are  parallel  to  x'  then  the  weight  of  that  edge  will 
equal  1  and  the  weight  of  the  other  edges  will  equal  zero. 


^x'  —  \^l,x'\^l—3,x'  T"  \@2,x'  1^2— 3, x' 


(38) 


where 


f*i— 3,^ 

h-3,x' 


n,x 

hx 


'  +  ^x' 
'  +  ^3,x' 


sign(/3ij;c/)Ft^_3^/ 1  Anm,l—3,x'  |  (dcfa— 3^x'  3,x',lim) 

sign(/32;ic/)Ft2_3?:c/ 1  A^m?2— 3,^' |  (dq.2— 3,x'  dc\2— 3,x',lim) 


dq  1-3, x' 
dq  2-3, x' 

dqi—3^x',lim 

dq_2— 3,x'  ,lim 

dqn:x' 


Ri-3,*'(qi  -qs) 

R2-3,^(q2  -  qs) 


<F  minmod 


$  minmod 


2dqi^/ , 
2dq.2,x/  5 


2dqi—3  x'i  2dq.3,a;/5  ^  (dqi^x'  H-  dq3x'^) 
2dq2-3,^5  2dq3jX/,  ^{dq2:X'  +  dq3 >a./) 


Ax  Hn,xrV qn:GGnx' 


(39) 

(40) 

(41) 

(42) 

(43) 

(44) 

(45) 


11  -38 


RTO-EN-AVT-186 


dMNATO 
WP  I  OTAN 


Challenges  to  Computational  Aerothermodynamic 
Simulation  and  Validation  for  Planetary  Entry  Vehicle  Analysis 


Note  again  that  $  is  the  same  auxiliary  limiter  used  in  Option  1. 
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